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Abstract 


Two sets of incremental objective stress measures and incremental non- lin ear strain mea- 
sures are compared for large deformation, dynamic, elasto-plastic s contact analysis with 
friction. These two sets are: (i) Jaumann stress rate tensor and incremental Green- 
Lagrange strain tensor, (ii) New incremental objective stress tensor and incremental log- 
arithmic strain tensor. 

A finite element computer code in FORTRAN is developed for the analysis of three 
dimensional, large deformation, dynamic, elasto-plastic, contact problem, including the 
friction at the contact interface. Elasto-plastic behaviour is modelled by an associated 
flow rule based on the von Mises yield criterion and power law type isotropic hardening. 
The material and geometric non-linearities are handled by using the updated Lagrangian 
method. A node-to-segment interface model is used for developing the contact stiffness 
matrix. Contact constraints are imposed using Lagrange multiplier method. Coulomb’s 
friction law is used for calculating the friction forces. A finite element - finite differ- 
ence scheme is used for spatial and temporal discretization respectively and Modified 
Newton-Raphson iteration method is used to obtain the solution of the non-linear govern- 
ing equation. To improve the convergence characteristics of Newton-Raphson iterations, 
under-relaxation technique and line search method are used. The finite-difference scheme 
used for time integration is Newmark’s algorithm. Unloading is incorporated to include 
the effects of global or local unloading during the analysis. Skyline scheme of assembly 
and static condensation scheme are used to reduce the computational resources required 
for the analysis. 
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Chapter 1 


Introduction 


Impact phenomenon occur in many engineering situations. The basic characteristics of an 
impact phenomenon are very brief duration, high force levels reached, rapid dissipation 
of energy and large accelerations and decelerations present. In many cases, the bodies 
undergo plastic deformations during impact. Sometimes frictional effects are also impor- 
tant during impact. A rigorous analysis of the impact problem requires the consideration 
of all these complex aspects. The analysis of the problems like the autmobile crashes 
involves the use of dynamics, elasto-plasticity and large deformation in a contact analysis 
scenario. 

Closed form mathematical solutions axe extremely difficult to develop even for simple 
impact problems. As the impact scenario becomes more complex, as in the case of elasto- 
plastic materials, arriving at a closed form solution is nearly impossible. This necessitates 
the use of approximate numerical tools for the analysis, of which finite element analysis 
is arguably the most versatile option. 

1.1 Review of Literature 

The classical incremental theory of elasto-plasticity [1] is generally taken to be the basis 
for predicting elasto-plastic response. The constitutive equation for these materials is 
derived by combining the stress-strain relations for elastic and plastic responses. The 
constitutive relation for a plastic response requires three fundamental characteristics of 
the material: a yield criterion, a flow rule and a strain hardening relationship. The yield 
criterion demarcates the elastic behaviour from the plastic behaviour. It specifies the 
state of multi-axial stress corresponding to the start of the plastic flow. Various yield 



criteria have been proposed by Tresca, von Mises etc. [2], where von Mises criterion is 
found to be the simplest and the most accurate in terms of the predicted values. The 
flow rule relates the plastic strain rate to the current stress. It is represented in terms 
of a plastic potential, analogous to the elastic potential, and is termed as the associated 
flow rule if the potential is taken as the yield function. The hardening rule specifies how 
the yield condition is modified during the plastic flow. For the purpose of analysis, the 
non-linear hardening relationship is usually modelled in terms of a power law. 

Unlike the case of small deformation of elastic materials, the governing equations 
of elasto-plastic materials, which undergo large deformation, contain two types of non- 
linearities: (i) geometric non-linearity associated with non-linear strain displacement re- 
lations and, (ii) material non-linearity in the form of non-linear stress-strain relations. 
For a contact problem, there is an additional non-linearity: Boundary non-linearity due 
to the non-linear boundary conditions at the contact interface. An incremental formula- 
tion like the Updated Lagrangian Formulation is the most convenient way to handle such 
non-linearities unlike the Total Lagrangian Formulation. The advantage of the former is 
that, if the increment size is sufficiently small, linearized incremental equations can be 
used. However, when the overall deformation is large, the number of increments becomes 
prohibitively large. In such a situation, one has to resort to one of the various itera- 
tive schemes available for solving the non-linear incremental equations. The prominent 
among such schemes is the Newton-Raphson method. In this method, the stiffness matrix 
is updated in each iteration of the increment. Thus, we work with the tangent modulus 
at that instant for each iteration. This method has two shortcomings. The first being 
that there is chance of divergence if the slope of force-displacement curve increases with 
the increment. Secondly, it takes more computational time as the stiffness matrix has to 
be computed afresh in each iteration [3]. Thus, the Modified Newton-Raphson method in 
which, the stiffness matrix is kept constant during the increment, is adopted. 

The dynamic analysis is very basic to studying impact problems. There are various 
ways of performing the time discretization in this category of problems of which, the 
direct integration methods are prominent. These methods include the central difference 
method, Houbolt method, Wilson-0 method and Newmark method [4]. The first method 
is conditionally stable one, while the rest of them are unconditionally stable. A study of 
the stability of Newmark algorithm in non-linear dynamics has been done by Hughes [5]. 
He has studied some potential problems where the stability of the algorithm fails. He 



has suggested that it can happen in a non-linear problem if the discrete internal energy is 
negative. Analysis of accuracy and speed of different algorithms is presented by Bathe [4]. 

It has been concluded that although central difference method is beneficial in terms of 
speed and cost, the trapezoidal method of Newmark’s algorithm is best in terms of stability 
and accuracy. 

1.1.1 Elasto-plastic Finite Element Analysis 

As stated earlier, the classical theory of elasto-plasticity [1] is generally taken to be the 
basis for predicting the elasto-plastic response. Argyris et al. [6] presented one of the 
earliest formulations and solution techniques for elasto-plastic problems. The effect of 
plasticity was simulated by “initial loads” and the matrix method was used to discretize 
the linearized governing equations. Two approaches, the direct incremental method and 
the iterative incremental method, were discussed for obtaining the solution. The early 
applications of the iterative incremental method [3] used the initial stiffness method. 
However, this approach is not suited when the overall deformation is very large as it 
results in very slow convergence, or in many cases, even divergence. 

Nagtegaal and Jong [7] dealt with the computational aspects of elastic-plastic large 
deformation analysis. A modified set of governing equations were formulated, which were 
applicable for small strain and large rotational increments. Special attention was paid 
to the integration of these equations in the deformation history. It is predicted that 
explicit methods, like the tangent modulus procedure, are conditionally stable whereas 
the implicit schemes, such as the mean normal method, are preferable. The main difficulty 
with the three dimensional analysis is the vast computational memory that it required 
for carrying out the analysis. But, with the advent of better machines and evolution of 
more efficient algorithms this problem has been solved and more and more such analyses 
are being acrried out. 

Since the elasto-plastic constitutive relation is generally given in the differential form, 
the constitutive function must be integrated to obtain the incremental stress. A major 
research effort has been devoted in recent years to the developement and critical assess- 
ments of stress updating or integration schemes for the rate (incremental) equation of 
plasticity. The procedures commonly used are the simple and robust Euler forward inte- 
gration technique, the implicit Euler backward integration technique and the radial return 



method. The trend has lately been shifting towards the use of consistent schemes which 
improve the convergence characteristics of the equilibrium iterations [8]. The recent works 
have clearly emphasized the use of implicit integration methods [9, 10, 11] in view of their 
superior stability and convergence properties [12]. A one step, fully implicit method of 
backward Euler integration scheme for large deformation of beams has been developed by 
Gendy and Saleeb [13]. 

Bathe [4] has briefly reviewed the different approaches used for obtaining the solution 
of general non-linear problems. Total and updated Lagrangian descriptions of motion 
have been used to formulate the incremental equations of motion. It has been noticed, for 
the problems solved, that the two formulations give identical results, provided appropriate 
constitutive relations are used. One of the areas of controversy in analyzing large defor- 
mation is the choice of appropriate objective stress measures. Various objective stress 
measures have been discussed by Bathe [4] and Chakrabarty [14]. These measures have 
been discussed in detail in section 1.1.4. The classical theory of elasto-plasticity assumes 
that the incremental strain tensor can be additively decomposed into the elastic and plas- 
tic parts. While this may be true for small incremental deformation, this assumption 
is no longer valid for large incremental strains and rotations, unless appropriate strain 
measures are used [15, 16] . 

1.1.2 Dynamic Finite Element Analysis 

Bathe et al [17] were probably the first to propose a large deformation dynamic finite 
element formulation. Total and updated Lagrangian descriptions of motion were used 
to formulate the incremental equations of motion. A detailed finite element analysis of 
non-linear static and dynamic response was also carried out by Mondkar and Powell [18]. 
The purpose of this paper is to review the theoretical and computational techniques for 
non-linear structural analysis, and their applications in a general computer program. A 
Lagrangian frame of reference is used and a variational form of the incremental equations 
of motion for large displacement, finite strain deformation is incorporated in the finite 
element program. It is found that no specific time integration scheme is optimal for all 
types of non-linear behaviour. It is inferred that for large load steps it is essential to use 
iterations to obtain accurate results. 

Direct integration methods are the most commonly used methods for solving the dis- 



cretized and linearized governing equations of motion. Dynamic equilibrium including the 
effect of damping forces is sought at discrete time instants instead of at all times. The 
advantage is obvious: all solution techniques employed in static analysis can probably be 
used most effectively in direct integration. The selection of any one direct integration 
method depends upon considerations of accuracy, stability and cost. The importance of 
using equilibrium iterations in dynamic non-linear analysis was first realized by Bathe [4]. 
An extensive discussion of these and other methods can be found in Bathe [4]. 

More recently, Gendy and Saleeb [19] found that, for a more comprehensive perspec- 
tive, work on the non-linear dynamic response of plates and shells calls for detailed studies 
of several important factors. These include the effect of large spatial rotations on the ge- 
ometric stiffness and inertia operators, accurate updating procedures for nodal rotations 
and associated angular velocities and accelerations, as well as material inelasticities (es- 
pecially for finite strains). Several of these issues were examined in conjunction with 
a newly developed mixed finite element formulation for plates and shells. To this end, 
and restricting the scope to the case of large overall motions but small strains, low-order 
displacement/strain interpolations are utilized, together with a radial return algorithm 
(backward Euler-integration scheme) for plasticity effects. The Newmark implicit scheme 
has been employed to integrate the semi-discrete equations of motion. 

1.1.3 Contact Analysis 

The analysis of contact has been a fascinating area since the initial days of the study 
of mechanics itself. A contact system may contain only one or more than one contact 
body, and contact-impact may occur between two bodies or within a single body. The 
behaviour of the contact system is governed by three main groups of equations, i.e the 
equation of motion, constitutive equations and initial and boundary conditions. The 
boundary conditions can be classified into two categories: prescribed boundary conditions 
and unkn own boundary conditions. The prescribed boundary conditions are deformation 
dependent. Among other possible deformation dependent boundary conditions, contact- 
impact conditions are an important class. Many investigators tried to determine the 
stresses and displacements in two elastic bodies coming in contact. In three dimensions, 
the problem of two contacting elastic bodies was first formulated and solved by Hertz [20] 
under several restrictive assumptions. Hertz assumed that the geometry of general curved 



surface in the vicinity of contact can be described by quadratic terms only. Neglecting 
friction and assuming that the bodies in contact deform as though they were elastic half- 
spaces, Hertz used the Boussinesq solution to deduce that the contact pressure distribution 
has to be elliptical to produce an elliptical area. For the contact between two spheres, 
the contact area is circular. Research work based on Hertz’s theory up to 1982 has 
been reviewed by Johnson [21]. Although Hertz evaluated the dimensions of contact 
area and stresses at the contact surface, he presented only a speculative sketch of the 
subsurface stresses. The stresses beneath the contact surface were analyzed later. For 
circular contact area of radius a, it is observed that the maximum shear stress occurs at 
a depth of 0.57a [21]. Later on Muskhelishvili [22] and Gladwell [23] solved some small 
deformation elastic contact problems using the analytical techniques. 

With the development of computers, people resorted to numerical methods. Finite 
element method has been the most widely used. In earlier studies, contacting bodies 
were assumed to consist of linearly elastic materials and hence were assumed to undergo 
only small deformation. The loading was assumed to be proportional. The node-to- 
node contact interface model was used for discretizing the contacting bodies and the 
contact conditions (constraints) involving nodal displacements and forces were applied 
by modifying the combined stiffness matrix or the flexibility matrix of the contacting 
bodies. The investigators who used this model are Ohte [24], Gaertner [25], Francavilla 
and Zienkiewicz [26], Sachdeva and Ramakrishnan [27] etc. 

Kalker and Randen [28], Hung and de Saxce [29] and Mahmoud et al. [30] used the 
technique of mathematical programmming for solving the frictionless small deformation 
elastic contact problem. They minimised the total strain energy to find the contact area 
and the contact pressure distribution. 

When the deformation is large, there is a possibility that a pair of nodes presently in 
contact, may undergo a large relative displacement. If a node-to-node interface model is 
used, re-meshing of the bodies becomes necessary for carrying out the subsequent analysis. 
To avoid frequent re-meshing, one must use the node-to-segment interface model, in which 
one node of a contacting body comes in contact with a segment of the other. When a 
node-to-segment interface model is used, the contact conditions acquire a complex form, 
when expressed in terms of nodal variables. To impose these conditions, an additional set 
of finite element equations, involving the contact stiffness matrix, needs to be developed 
using the principle of virtual work of the contact forces. The methods, which are normally 



used for implementing the contact constraints are Lagrange multiplier method, penalty 
funtion method and transformation matrix method. Bohm [31] presented a comparison of 
different contact algorithms with applications. A survey of contact algorithms is given by 
Bourago [32]. 

Kikuchi and Oden [33] have provided the mathematical foundation for the Lagrange 
multiplier method. The work of Bathe and Chaudhary [34] seems to be the first at- 
tempt to develop an expression for contact stiffness matrix using this method. This was 
done for planar and axisymmetric problems assuming the contact segment to be straight. 
Lagrange multiplier was used to impose the contact constraints including friction forces 
based on Coulomb ’s friction law. Non-linear contact kinematics was used for developing 
a consistant contact matrix for two dimensional problems by Wriggers and Simo [35] and 
by Parisch [36]. They used both the Lagrange multiplier method and penalty function 
method. The transformation matrix method for two dimensional problems was proposed 
by Chen and Yeh [37]. A formulation for non-linear frictional model was proposed by 
Gallego and Anza [38] using perturbed Lagrangian functional. All the above formulations 
were used for static large deformation elastic contact problems. Some of these formula- 
tions were extended to elasto-dynamic contact problems [4, 35] and dynamic elasto-plastic 
problems [39]. A literature survey of contact dynamics modelling is presented by Gilardi 
and Sharf [40]. Sung and Kwak [41] and Bittencourt and Creus [42] have presented the 
analysis of dynamic, large deformation contact problem including the effects of friction. 
Master’s thesis by Werner [43] is an attempt to compare the experimental results and 
the results predicted by a commercial finite element software for an elasto-plastic impact 
problem. 

Contact-impact algorithms, which are sometimes called as slideline algorithms, are the 
computationally time-consuming parts of many explicit simulations of non-linear prob- 
lems, becuase they involve so many branches, and thus are not amenable to vectorization, 
which is essential for speed on supercomputers. Belytschko and Neal [44] introduced the 
pinball algorithm, which is a simplified slideline algorithm and is readily vectorized. Its 
major idea is to embed pinballs in surface elements and to enforce the impenetrability 
condition only to these pinballs. It can be implemented either by a Lagrange multiplier 
or penalty method. Malone and Johnson [45] developed a parallel contact algorithm and 
implemented it in a non-linear dynamic explicit finite element program to analyze the 
three-dimensional shell structures. The contact algorithm accounted for initial contact, 



sliding and release through the use of parametric representation of motion of points lo- 
cated on the surface of structure combined with a contact surface representation, which 
approximates the actual surface by means of triangular search planes. Brown et al. [46] de- 
scribed a general strategy for parallelizing solid mechanics simulations. Such simulations 
often have several computationally intensive parts, including finite element integration, 
detection of material contacts and particle interactions if smoothed particle hydrodynam- 
ics is used to model highly deforming materials. Their strategy is to load-balance each 
of the significant computations independently with whatever balancing technique is most 
appropriate. The main advantage is that each computation can be scalably parallelized. 
The drawback is the data exchange between the processors and the extra coding that 
must be written to maintain multiple decompositions in a code. 

The problems of undesirable oscillatory solutions and stability of time-integration 
schemes have been tackled by quite a few people. Taylor and Papadopoulos [47] address 
the problem of oscillatory solution along the contact interface. The standard second order 
implicit integrators of the Newmark family, used in earlier works, were investigated by 
them. In their work, control of these oscillations is attempted by introducing an artificial 
bulk viscosity for the compressive waves in each body. Laursen and Chawla [48] proposed 
an energy conserving algorithm to control the oscillations for the frictionless dynamic 
contact problems. It is shown that a Lagrange multiplier enforcement of an appropriate 
contact rate constraint produces these conservation properties. In particular, the ability of 
the formulation to produce accurate results, where more conventional integration schemes 
fail, is emphasized by numerical simulations. Armero and Petocz [49] presented a new 
dissipative time-stepping algorithm to improve the stablity of time integration scheme. 
They showed that although the traditional time-integration schemes were stable for most 
of the linear problems, instability would creep in as soon as we try to use them in non- 
linear cases. Given these considerations, they developed implicit time-stepping algorithms 
for contact problems that posess unconditional (energy) stability in time and lead to stable 
enforcement of contact constraints. 

Quite a few researchers have proposed new techniques. Ayari and Saouma [50] devel- 
oped a new formulation for contact problems using the concept of fictitious forces. Con- 
trary to the most existing models, this one involves very few matrix decompositions and 
thus is computationally inexpensive. The model was applied to fracture mechanics based 
analysis of a cracked dam under seismic excitation. They also showed that their model 



could be applied to a vast majority of contact problems and was not only restricted to 
fracture phenomenon. Simo and Laursen [51] extended the augmented Lagrangian fomula- 
tion, which has certain advantages over traditional Lagrangian formulation, to problems 
involving frictional contact. Kane et al [52] proposed a nonsmooth contact algorithm, 
which is capable of dealing with multibody nonsmooth contact geometries for which nei- 
ther the normals nor the gap functions can be defined. Farahani et al [53] presented a 
solution method for contact/impact problems, which is different from Lagrange multiplier 
and penalty methods. This method is based on the stiffness matrix transformation and 
eliminating the normal degree of freedom of contactor node. The method is absolutely 
general and can be used in static as well as dynamic non-linear problems. In this case, 
the number of unknowns does not increase and the solution scheme needs no modification 
as is the case with prior methods. Besides, the boundary conditions are satisfied exactly. 

1.1.4 Objective Stress Measures 

In large deformation elasto-plastic problem, one comes across large rotations also besides 
large displacements and strains. In updated Lagrangian formulation, this creates difficulty 
in updating of stress components. Note that the constitutive equation gives the stress 
increment only due to incremental straining of the material. If the Cauchy stress tensor 
is used as the stress measure, then one cannot simply add the incremental stress to the 
Cauchy stress components at time t to obtain their values at time t-f At. This is because, 
the change in Cauchy stress components also occurs due to the incremental rotation. Thus, 
the Cauchy stress tensor is not objective in the sense that its components do change even 
when the increment consists of pure rigid rotation. The stress measures whose increment 
(or rate) depends only on the incremental deformation (or rate of deformation) are called 
as objective stress (or stress rate) measures. 

The most commonly used objective stress measure in the literature as well as in finite 
element software on large deformation elasto-plastic problems is the Jaumann stress rate 
measure. It is obtained by eliminating the rate of change of Cauchy stress component 
due to rotation from the total rate. But it uses only the infinitesimal rotation tensor 
to represent the incremental rotation. Therefore, when the incremental rotation is large, 
it is not expected to give accurate results. Further, it is well known that the use of 
Jaumann stress rate measure predicts an oscillatory solution when the incremental shear 



is large [54]. Other objective stress rate measures like the Green-Naghdi rate and E-rate 
have been proposed ([54] [55]). But they are difficult to use. A simple stress updating 
procedure has been proposed by Varadhan [56] which accounts for the change in Cauchy 
stress components due to incremental rotation. In this approach, the incremental stress is 
evaluated in a frame which rotates with the particle rather than in a fixed frame. While 
transforming the stress components from the material frame to the fixed at time t + A t, 
the transformation due to incremental rotation is taken into account. This incremental 
stress measure has been called New incremental objective stress measure in this work. 
This stress measure, along with other objective stress measures, has been discussed in 
detail in section 2.3. 

1.2 Objective of the Present Work 

A finite element code for the analysis of 3-D, large deformation, elasto-plastic, dynamic, 
contact problem with friction has been developed by Jayadeep [57]. It uses the Jaumann 
stress rate measure and incremental Green-Lagrange strain measure. The objective of 
the present work is to incorporate the new incremental objective stress measure and 
incremental logaritmic strain measure in this code and to carry out the comparitive study 
of these two sets of stress and strain measures. 

The salient features of the code are as follows. Since we have large deformations, 
including slipping at the contact interface, a node-to-segment contact model is used in 
contact analysis. A method outlined by Zhong [39] is used for calculating the contact 
stiffness matrix. Lagrange multiplier method is used for applying the contact constraints. 
Coulomb’s friction law is used for modelling the frictional effects at the contact interface. 
Elastic behaviour is modelled using generalised Hooke’s law and the plastic behaviour is 
modelled by an associated flow rule based on von Mises yield criterion and power-law 
type isotropic hardening. The total strain is obtained by adding the elastic and plastic 
components. Incremental analysis is carried out using the updated Lagrangian method [4]. 
Modified Newton-Raphson iterative method is used for solving the resulting finite element 
equation. In each Newton-Raphson iteration, contact iterations are carried out to find 
the contact reactions. A finite difference scheme is used for carrying out the discretization 
in time. Newmark’s algorithm, which is found to be the best in terms of stability and 
accuracy among the various finite difference schemes, is used for this purpose. Unloading 



scheme is incorporated to take care of the local unloading during the analysis. To reduce 
the computational resources required for performing the analysis, static condensation of 
stiffness matrix for contact and non-contact nodes and skyline scheme of assembly are 
used. Divergence handling techniques, viz., under relaxation and line search, are used for 
improving the convergence characteristics of Newton-Raphson iterations. 

1.3 Structure of the Thesis 

Chapter 2 deals with the mathematical formulation of the governing equations for dy- 
namic, large deformation, elasto-plastic problem. The development of the finite element- 
finite difference scheme based on these governing equations and the numerical scheme 
used are discussed in chapter 3. Chapter 4 presents the contact formulation, where in 
the formulation of the contact stiffness matrix and the calculation of contact reactions 
are discussed. Chapter 5 is devoted to the results of the comparison between the two 
sets of objective stress measures for elastic, elasto-plastic, single body as well as contact 
problems. The last chapter summarises the work carried out and the scope for future 
developments. 



Chapter 2 


Mathematical Modeling of Dynamic 
Elasto-Plastic Problem 


In a dynamic elasto-plastic problem, the deformations and strains in the body no longer 
remain small and hence the response of the body is not proportional to the load. The 
response of the body becomes non-linear. The non-linearities can be broadly classified 
into three: Geometric non-linearity, material non-linearity and boundary non-linearity. 
First two of these non-linearities occur in the analysis of large deformation of elasto-plastic 
materials. The most convenient formulation for such analysis is the Updated Lagrangian 
Formulation with the constitutive relations represented in an incremental manner. The 
third non-linearity occurs in the analysis of contact-impact problems and this is attributed 
to the non-linear boundary conditions at the contact surfaces. 

2.1 Introduction 

In the study of the deformation of a body subjected to external loading, often the original, 
undeformed and unstressed state of the body is used for the formulation of its governing 
equation. This method is known as the Lagrangian formulation. This formulation is 
convenient for small deformation problems, where the deformed configuration does not 
deviate from the original configuration significantly. Hence, the deformations are pro- 
portional to the loading and can be described by an infinitesimal or linear strain tensor, 
for which the strain - displacement relations are linear. On the other hand, for large 
deformation problems, one has to use a finite strain measure, which is expressed by a 
non-linear strain displacement relation. Furthermore, the equations of motion, expressed 



in the original configuration, depend on the deformation. This makes the governing equa- 
tions highly complex and too difficult to solve. In such cases, it is convenient to solve 
the problem in an incremental manner known as the Updated Lagrangian Formulation. 
In this formulation, the state of the body is updated in increments using the calculated 
incremental deformations and stresses. The response during At is calculated based on 
the configuration at time t, where At represents the increment in time, called time step. 
Hence, all the relations are to be in an incremental form rather than the total form. This 
methodology is particularly useful for elasto-plastic materials, since the constitutive rela- 
tions for these materials depend on the state of the body and therefore they are expressed 
in an incremental form. 

2.2 Kinematics of Finite Deformation 

The relative deformation gradient tensor at time t 4- At is defined by: 

t+A t% = = t+A txi,j (2.1) 

where t+At X{ and t Xj denote the position vectors of the particle at time t and t + A t 
respectively. When t+At F is not singular, the polar decomposition theorem ([58] , [59]) 
allows a decomposition of the form: 

t+Mpi. = t+At^ t+ At^ (2.2) 

where £+A£ R;k is an orthogonal tensor representing the material rotation and t+A tUkj is a 
positive definite symmetric tensor called the right stretch tensor. The right stretch tensor 
can be diagonalized by the following transformation to obtain the principle stretches 

t+At\ . 
tAj. 

t+A t = t+A lQik t+A tUki t+A tQji (2.3) 

where t+A t[Q] is orthogonal matrix and 

t+A tXi - t+A tUfi {no sum ) (2.4) 

The tensor t+A /R and the transformation matrix, t+A t[Q] are used in the development of 
the stress measure while the incremental logarithmic strain tensor is obtained from the 



2.3 Stress Measures 


The large deformation problems are normally associated with rigid body rotations. Hence 
the stress measures used in these problems should be objective, which means that they 
should vanish if the increment is a pure rigid body rotation. The most commonly used 
stress measure known as Cauchy stress tensor is not objective and hence it can not be 
used in a constitutive equation. There are numerous objective stress measures, of which 
some are discussed below. 

One of the commonly used objective stress measures is the Second Piola-Kirchoff 
stress tensor t+A *Sij which is related to Cauchy stress tensor t+At cr mn using the concept 
of equivalent work between the two configurations [4]: 


t+Atn 

t&ij — TTTT 


t+ Ats 


i t+At_ 
t+At ^ijTn 


ran t+At % j 7 n 


(2.5) 


Here t p and t+At p represent the densities at time t and time t + At respectively, t +At %i,m 
is the derivative of the position vector at time t, viz., t x i with respect to t+Ai x m . 

o 

Another frequently used objective stress measure is Jaumann’s stress rate tensor ffy. 
It is related to the Cauchy rate by the relation [4]: 


t a ij At = dij At - bikCfyk At) - bjkC&ik At) (2.6) 

^Cljk At = — (fAuij — tAuj,i) (2-7) 

Here, At represents the incremental spin tensor and t Auij denotes the derivative of 
the incremental displacement vector t Aui with respect to hj. 

Dienes [54] proposed a stress rate d based upon the rotation part of the deformation 
gradient: 

Cij — dij Qik&kj + Cikj&ik ( 2 - 8 ) 

where, 

Qij = Rik Rjk (2-9) 

is the angular velocity. The stress rate d is sometimes refered to as the Green-Naghdi 
rate. 

As pointed out by Metzer and Dubey [55], it is important that the stress measure 
used be compatible with the constituve equation in addtion to being objective. A new 
incremental objective stress measure (as against objective stress rate measures) to be used 
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Note: 

• The Fixed Frame is denoted by subscipt F. 

• The Material Frame is denoted by subscript M. 

Figure 2.1: Fixed and material frames 

in the generalized Hooke’s Law (2.16) and in the elasto-plastic constitutive equation (2.39) 
is developed below. 

Two Cartesian referance frames (see Figure 2.1) are used: 

1. The fixed frame (denoted by subscript ,F). 

2. A material frame (denoted by subscript ,M) which rotates and translates along with 
a material particle. 

The material frame is defined so as to coincide with the principal axes of the right stretch 
tensor at time t so that the initially orthogonal axes do not get skewed at time t + At. 

The first step in the derivation is the transformation of the components of the Cauchy 
stress tensor at time t from the fixed frame to material frame: 

‘cr%= ,+A tQ a ‘au ,+i ?<3ji (2-10) 

The increment in the components of the Cauchy stress tensor with respect to the material 
axes, tAtfij, is added to t<?ij to obtain the stress components at time t + At with respect 
to the material frame of reference: 

t+At a*f = *o}f +* A eg (2.11) 

The final transformation from the material frame at time t + At to the fixed frame gives 
the components of the Cauchy stress tensor at time t + At: 

= M tQu ,+A 1 ‘H i .‘ +A V," t+A lKm ,+ “Q„j 


( 2 . 12 ) 



where t+A tH has components with respect to fixed frame. 

The proposed stress and logarithmic strain measures satisfy the req uir ements of ob- 
jectivity and lead to a physically consistent application of the usual constitutive equation. 

2.4 Strain Measures 

The linear or infinitesimal strain tensor, normally used in the small deformation problems, 
is not a valid choice for the large deformation, elasto-plastic analysis. Green- Lagrange 
strain tensor is one of the most commonly used strain measure for large deformation 
analysis. 

The incremental Green-Lagrange strain tensor is a non-linear function of the displace- 
ment [4]: 

tAejj = “(tArtij -1- tAu^, -1- tAu^j tAu^j) (2.13) 

When the incremental deformation is small, the incremental linear strain tensor 

t Asjj = ~ (t Ai£jj + {A Uj t {) (2.14) 

can be used as the measure of deformation. 

These two strain measures occur in the virtual work expression at t + At and its 
transformation to time t respectively. The components of the Green Lagrange strain 
tensor are invariant under rigid body rotation of the material unlike the small strain 
components. The following are the disadvantages one would face using one of the above 
strain measures in a constitutive law: 

1. The solution obtained is dependent upon the size of the increment in the updated 
Lagrangian formulation unless the increment size is sufficiently small. 

2. The components of the strain tensors do not tend to infinte values when the prin- 
cipal stretches tend to zero. Therefore a constitutive law which ensures that the 
appropriate Cauchy stress components tend to negative infinity (as is physically 
realistic), even though the strain components remain finite, should be used. This 
difficulty can be avoided by using a strain measure whose components become minus 
infinity when the principal stretches become zero. 

The Logarithmic strain measure introduced by Dienes [54] is free of the above disad- 
vantages. The components of the incremental logarithmic strain tensor, in the material 



frame, are defined by: 


t Aeg = ln( t+ ^Xi) (no sum over i) (2.15) 

where the t+A f A* are the principal stretches defined by equation (2.4). The components 
with respect to any other frame can be obtained by the usual tr ans formation law. The 
logarithmic strain has the following additional advantage in elasto-plastic analysis. A 
loading test involving elasto-plastic deformation followed by elastic unloading reveals that 
the slope of the elastic unloading line is the same as that of the initial elastic loading line 
only when the true stress and the logarithmic strain measures are used in a constitutive 
law. 


2.5 Elastic Constitutive Equation 


The elastic constitutive equation used is the generalized Hooke’s Law relating the incre- 
mental objective stress tensor and the elastic part [Ae ei ] of the incremental logarithmic 
strain tensor: 


= c§ t , ,a 4 f 

(2.16) 

The tensor C? kl for isotropic case is given by: 


Cfj k i = A 6ij 5ki + 2 n SikSji 

(2.17) 


where A and // are Lame’s constants. 

The stress and strain measures used in a constitutive equation need not necessarily 
be energy conjugate with each other. However if it is so, the predicted response in a 
predictor-corrector scheme will be closer to the actual response. 


2.6 Theory of Plasticity 

An Elasto-plastic material model is characterized by the following: 

1. Yield Criterion. 

2. Strain hardening postulate. 

3. Flow rule. 

Baushchinger effect and the hysterisis loop are disregarded and the material is assumed 
to be isotropic in this work. 



2.6.1 Yield Criterion 


A law defining the limit of elastic behaviour under any possible combination of stress is 
called yield criterion. This law applies both to loading directly from the annealed state 
and to reloading of an element unloaded from a previous plastic state. The von Mises yield 
criterion which has been experimentally verified to predic the yield locus quite accurately 
is used. Yielding occurs when the equivalent stress t a eq reaches the value of the yield 
stress in uniaxial tension t a y . 

t or e g(ai j )= t a y ( 2 . 18 ) 

where, 



and t a' i j is the deviatoric part of t a i j. • 

2.6.2 Strain Hardening Postulate 

The size and the shape of the yield surface in a strain hardening material depends upon 
the complete history of plastic deformation since the prevoius annealing. A convinient 
mathematical formulation for strain hardening is obtained by assuming that the yield sur- 
face expands uniformly without changing its shape. This is known as isotropic hardening 
and the relationship may be stated as 

t<T y = fC< g ) (2-20) 

where, 

L (2 - 21) 

is the equivalent plastic strain determined from the path-line integration of an invariant 
of t de P if (plastic part of the incremental logaritmic strain tensor) and the function / is 
determined by the relationship between the true stress and the plastic strain in uniaxial 
tension/compression. In this work, the following power law is chosen to represent /: 

\ - \ = K(%r (2.22) 

where °a y denotes the initial yield stress of the material. 



2.6.3 Plastic Flow Rule 


The rule of plastic flow relates the current stress with the plastic strain increment. The 
plastic potential is a function defining the ratios of the components of the incremental 
plastic strain. The associated flow rule which assumes that the yield function and the 
plastic potential are the same is used in the present work. It states that 

ftp 

tdef = dX-^- . (2.23) 

where dX is a scalar to be found from the hardening relation. 


2.7 Elasto-plastic Constitutive Equation 

An elasto-plastic relationship between the incremental stress and strain tensors based on 
the von Mises yield criterion and isotropic hardening (as discussed above) is developed 
below. 

For an isotropically hardening material, the plastic potential is given by [60] 1 : 

F{ cr ij) F Zeq) = Peqi&ij) ~ °y( Pe eg) (2.24) 


where, 

F = 0 


(2.25) 


represents the yield criterion. The plastic potential F depends on the Cauchy stress tensor 
through the equivalent stress, and on the variable yield stress through the hardening 
parameter p e eq . 

The plastic part of incremental logarithmic strain tensor d<%f is obtained from the 
plastic potential using the flow rule (equation 2.23). Differentiation of equation (2.24) 
with respect to cry gives: 


dF _ 3 , 

doij ~ 2a ^ 


(2.26) 


Therefore, dX is determined as: 

i\ = *5, (2.27) 


Further, the hardening relationship and the yield condition can be used to express dX as: 


day da ^ 

dA= ir = ir 


(2.28) 


1 In this section, the superscript/subscript denoting time have been omitted for convenience. 



where, 


da, 




eg 


(2.29) 


is the slope of the hardening curve. Substitution of equations (2.26) and (2.28) in equa- 
tion (2.23) leads to the following constitutive equation: 




3 da. 


v 


—a 1 -- 

2Ha eq a * 


(2.30) 


This constitutive relationship between the deviatoric stress tensor and the plastic 
part of the incremental logarithmic strain tensor is not really convenient for the Updated 
Lagrangian formulation for which an incremental stress-strain relationship is needed. This 
can be obtained from equation (2.30) as follows: 


_ 3 a ij 9a eg 

< - w^5^ dx!u 

Note that, from equations (2.24) and (2.26), we get: 


(2.31) 


da. 


eq 


dF 


~ a ki 


da u da k i 2a eq ~* L ^ 2 ' 32 ^ 

Substitution of equation (2.32) in equation (2.31) leads to the following incremental plastic 
stress-strain relationship: 

(2-33) 

The incremental elastic stress-strain relationship (equations (2.16) and (2.17) in the 
inverted form) can be expressed as: 

de\f = ^[-v da kk 5 {j + (1 + v) da y ] (2.34) 

where, de ff is the elastic part of defj, E is the Young’s modulus and v is the Poisson’s 
ratio. Adding the two relationships, we get: 


deh = deff + d<$- 


'-ij 


v 


v r r , 1+ v s s , 9a J ij a l kl 
—6 {j S u + 15 j|- 


da k i 


(2.35) 


This is the incremental elasto-plastic stress-strain relationship needed in the updated 
Lagrangian formulation. However, it is the following inverse relationship which is more 
useful ([60] and [61]): 

daij = Cg ^4) (2.36) 

where, 

/ U i < /t' . rr'. . \ 

(2.37) 


„ pp n (, . , v . . 9 POjjQu 

Cijki ~ Wji + 1 _ 2u 5ij6kl 2(3 fi + H)a% 



Here, /z is one of the Lame’s constants (also called as shear modulus) and v is Poisson’s 
ratio. The Poisson’s ratio is related to the Lame’s constants A and fj, by the following 
relation: 

A 


v = 


2 ( a + ,y < 2 - 38 > 

Note that the stress increment appearing in equation (2.36) must be an objective 
stress increment in the sense that dcij must be a zero tensor if the increment is a pure 
rotation. The incremental stress measure used in the present work is the New Objective 
stress measure discussed in section 2.3. 

The relationship (2.36) has been derived assuming the increment size to be small. 
When a large size increment is to be used, the relationship (2.36) takes the following 
form: 

-t+ At 

‘CFMM,) (2-39) 

where, 

\ 


•C t %=2u(6 ik 6 il+ ^JL^SqSu- 


2(3 I/.+ (2 ' 40) 

Here, H has been replaced by the expression (2.29) and the left subscript t has been added 
to make it explicit that these quantities are to be evaluated at time t. 


2.8 Incremental Updated Lagrangian Formulation 

The objective of the updated Lagrangian formulation is to establish dynamic equilibrium 
in the configuration at time t + At when all the static and kinematic variables at time t 
are known. The principle of virtual work requires that [4]: 

J t+Mp t+Mjli 5 (‘ +A ^) dt+AtV + [ t+Atv t+Ata v 6 ( t+At£ a) dt+Aty = t+MR ( 2 - 41 ) 

where, 

t+At p = density at time t + At 

t+At Ui = acceleration at time t + At 

S( t+Ai Ui) = virtual displacement at time i + At 

t+At V = volume at time t + At 

t+At &ij = Cauchy stress tensor at t + At. 

The virtual linear strain tensor 8( t+At Sij) is defined as: 


^ lj) 2 \d t+At Xj d t+At Xi J 


(2.42) 



Further, 


,+AtR = h* v ,+a%i 5 ( t+A ^o <* HA v + / ,+A, i:t i ) s i(‘ +a v j )d ,+i, s+ 

/ +Atfc 1+ ^)oJ(‘ + ^)^ (2.43) 

where, 

t+A %i = body force per unit volume at time t + At 
t+A \ti)s = applied traction per unit area at time t + At 
t+At Sp = surface at time t + At with traction specified 
t+A %ti)c = contact traction per unit area at time t + At 
t+At Sc = contact surface at time t + At 

The main difficulty in applying equation (2.41) is that the configuration at time t+ At 
is unknown. Therefore, the virtual work expression at time t + At is transformed to an 
integral over the volume at time t. It is assumed that the external load term (2.43) is 
deformation independent for the formulation of governing equation. The expression (2.41) 
after the transformation becomes [4] : 


J tv V t+A % 5( t+A< Ui) d l V+ J t+A £% <K t+A teii) d*V = t+At R 

where the virtual Green-Lagrange strain tensor d( t+A *etj) is defined by: 

* ( t + At x _ 1 \d5P*ui ' 

1 t6lj) 2 d t+At Xj d t+A % \d t+At Xi d t+At Xj ) 

The second Piola-Kirchoff stress tensor can be decomposed as: 


(2.44) 


(2.45) 


t+A &; = ISij + t ASij = % + t ASij (2.46) 

Since, 

5( ,+At u i) = ifu j + ,Aui) = S( t An,) (2.47) 

the virtual Green-Lagrange strain tensor can be decomposed as: 

S( t+A t e ij) = = KAsij + tAriij) (2.48) 


tAsij — g {tAuij + tA Uj i i) (2.49) 

tAtiij — —{tAiikj fAukj) (2.50) 


where, 



Therefore, equation (2.44) can be written with incremental decomposition as: 


J ty V t+A % 5(Aui) d*V + t ASij 5( t AEij) dW + 

J tv tAS'y SitAriij) d t V+ t a ij S^A^) d*V + 

J % 5( t Aeij) d l V = t+At R (2.51) 

The above equation is simplified by neglecting the third integral, which is a higher order 
term, and approximating t ASy as t C^ l 

j tv t P t+A %5( t Au i )i‘V + j iv t C^ t Ae kl S( t Ae il )d t V + 

J V 0 JfiAijj,) d‘V + J 'ct,, 4( t A £li ) d‘V = ,+a *R (2.52) 

The linearized equation, when solved, will yield only approximate displacement, ve- 
locity, acceleration, strain and stress fields. The approximate quantities are denoted by 
a right superscript (1). The error due to the approximation involved is calculated from 
equation (2.41) as: 

Error = ,+At R- [ ,+A Vg> i( ,+A ^g>) d‘+ A V< l > - 

Ji+AVO) 1 3 

l A V (1) t+Atp{l) t+A ^ X) ^( t+At 4 1) ) d t+A V (1) (2.53) 

This error is generally minimised by a predictor - corrector scheme. 



Chapter 3 


Finite Element-Finite Difference 
Formulation 


The governing equation (2.52) contains volume and surface integrals, which can not be 
solved analytically for most of the physical problems. In many of the cases of interest, the 
volume and surface may not be expressed mathematically. This necessitates the use of 
numerical methods for the solution of the governing equation. Since it approximates the 
field variable and geometry simultaneously, the Finite Element Method is probably the 
most powerful method for solving the complex boundary value problems. The integration 
in the time domain, which is an initial value problem, is carried out efficiently using 
a Finite Difference Scheme. This chapter presents the development of a finite element 
- finite difference scheme for solving the linearized governing equation of the dynamic 
elasto-plastic problem. A numerical iterative technique for minimizing the error term 
given by equation (2.53) is also presented. 

3.1 Matrix Notation 

Matrix notation is used in the development of finite element — finite difference equations 
for the convenience in computer implementation. 

The components of the strain tensors *A % and t A %• are represented in the array 
form as follows: 


== {t Ae XI , tAe^y, t Ae zz , 2^Ae x ^, 2tA€y Z , 2 t Ae xz ) 
t{A 77}^ = {t Au^, t tAu jZ , jAi^a;, . . .} 


(3.1) 

(3.2) 



} {t^ e xx> j/t/i t^ e z«) tAe£ y , tAe^ 2 , tAe£ 2 } (3.3) 

Note that the shear strain components in the array t {Ae} 7 contain a factor of 2, which is 
a convention followed in most finite element formulations. The components of the Cauchy 
stress tensor Vy are written as: 

{ t } { &yyi ®zzi &xy> ®yzi &zx} ( 3 - 4 ) 

Two matrix forms of the tensor t C EP l are needed: 

1. t [C EP ] is used in the constitutive relationship t {A5'} = l [C EP ] t{Ae} 

2. t [C Epl ] is used in the constitutive relationship t {Aa M } = f {C EP '\ t{Ae L } 

These matrix form are given by: 


U/^EPl 


t\r<EP' 


In the above equations, 


: |C E ] <{o} ‘{a} T [C E ] \ 

1 *»+ ‘{a } T [C E \‘{a}) 

(3.S) 

J *H+ ‘{a} T [C E> ] *{a}) 

(3.6) 

— ^ t r n 

2‘cr„ * 

(3.7) 


(3.8) 


The array t {o'} T represents the deviatoric part of the Cauchy stress at time t. For an 
isotropic material, the matrices [C E> ] and [C E ] for the 3-D case are given by: 


[C B ] = 


(1+ u){l- 2v) o 


(1 + i/)(l - 2v) 
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(3.10) 



Equation (2.52) can be written in the following form owing to the symmetries of 


t Aeij, tArjij and t a ij : 


J v ^.{Auf) V ,+A *{«} i'V +.J iv J(«{ A<f) *[C EP ] ,{Ae} d‘V 

+ f, v S('l^} T )‘{T\ t {&n}d‘V + J' V 6( t {Ae} T )‘{o}d‘V = t+i *R( 3.11) 


The matrix [T] is given by: 


where, 


*[E] 0 0 

*[T\ = 0 *[E] 0 

0 0 *[E] 


O xx 0 X y O' xz 

[E] — t O X y t Oyy t Oy Z 

*0 XZ <Jy Z <J ZZ 


(3.12) 


(3.13) 


3.2 Finite Element Formulation 


The domain is discretized into a number of elements and the incremental displacement 
field is approximated over each element by: 

f > 

t Au 

t{Au} = i t Av = ‘[$3 t{Au} e (3.14) 

t Aw 

t j 

Here, the element incremental displacement vector is given by: 

({Ait} er = {tAu 1 , [A?/, t Aw l , (A« 2 , . . .} (3.15) 

where the quantities ( A u\ tAv', t Aw' stand for the unknown incremental displacements 
of node i in x, y and z directions respectively. The matrix *[$] is defined by: 

*[«] = *{®2} r ( 3 ' 16 ) 

where, 

‘{®i} T = pV 1 ,0,0,‘Af a> 0 > 0,...} 

‘{®2} t = {0,*1V 1 ,0,0,*« 2 ,0,...} 

‘{«s} T =• {0,0,'JVr, 0,0, %,...} 


(3.17) 



The W;, which are functions of (x,y,z), are called as shape functions. In the present 
work, 8-noded brick element with trilinear shape functions is used. The functions Wj do 
not depend explicitly on t. However, the left superscript t is used to em ph asiz e the fact 
that the shape of element (nodal coordinates) changes with tim e. 

The acceleration can similarly be expressed as: 


t+At ii 




t+At v > 


t+A1 u 


4 [0] t+At {u} e 


(3.18) 


The strain field is expressed in terms of the nodal displacements by differentiating 
(3.14) and using the expressions (2.49) and (2.50) as: 


t {Ae} = WtiAuY (3.19) 

.{At;} = l [B„] t {A u }* (3.20) 


where, 


and 


w 


WI 

w; 

w; + 
[w;+w; 


W T = {*{$1},*, W*.-} 


(3.21) 


(3.22) 


Using equations (3.14), (3.18), (3.19) and (3.20), the contribution to the integral (3.11) 
over a typical element with volume V e is: 


5(,{Ati}' r ) (J *[*F V *[*] i’ v “) ,+a '{fi }' + 
i( t {A«}' T ) (J^ ‘[BiF '[C BjP ] *[Bl] d*v) .{Au}' + 
i( t {A«}' r ) (/ v> ‘[B„] T ‘[T] ‘[B w ] dV) ,{A«}'+ 

< S(,{A»}‘ r ) (/ y< ‘[B t ] T *W W) = i(«{A«}'") (3-33) 

The contribution to the term t+At R is expressed in terms of the elemental external force 
vector t+At {F} e using a standard procedure [62]. 



Since the variation in the displacement vector is arbitrary, the above equation can be 
written as: 

W t+At {u} e + \KY *{A«}« + ‘{/} e = t+At {F} e (3.24) 

where the element mass matrix t [M} e is defined as: 

*[M] e = f ‘[SfVWV 8 (3.25) 

and the element stiffness matrix *[. K ] e is expressed as: 

W = W+ ‘[W (3.26) 

where, 

^KlY = [ (3.27) 

jtyc 

‘[■Katl] 6 = [ (3.28) 

Jtye 

The element internal force vector is: 

V} e = f tve t {B L ] Tt {a}d t V e (3.29) 

The element mass ( f [M] e ) and stiffiiess ( 4 [A] e ) matrices along with the element force 
vectors t {f} e and t+At {F} e are assembled to obtain the global equation: 

t [M] t+At {U} + *[K] t {AC/} + *{/} = t+At {F} (3.30) 

where, 

t [M] = global mass matrix, 
t [K] = global stiffness matrix, 

*{/} = global internal force vector at time t, 
t+At {F] = global external force vector at t + At . 

The equation (3.30) represents a system of ordinary coupled second order linear dif- 
ferential equations. The next section describes the technique to convert it into a system 
of algebraic equations. 

3.3 Finite Difference Scheme 


The number of efficient methods for the solution of a large system of coupled ordinary 
differential equations such as those arising from a finite element formulation is limited. A 



popular technique is the direct integration of equation (3.30), where the dynamic equilib- 
rium is sought only at discrete time intervals instead of at all times t. The assumptions 
on the variations of displacement, velocity and acceleration within each time interval 
determine the accuracy, stability and cost of the solution procedure. 

Two broad classifications exist in the direct integration methods: Explicit and Implicit 
integration methods. Explicit integration methods are based upon the equilibrium con- 
ditions at time t for obtaining the solution at time t + At, resulting in such procedures 
being only conditionally stable. One main advantage of explicit schemes is that factor- 
ization of the effective stiffness matrix is not necessary. Implicit integration methods use 
the dynamic equilibrium conditions at both time t and t + At for obtaining the solution 
at time t + At. Since the implicit time integration methods are unconditionally stable, 
the time step used can be several orders of magnitudes larger, but the factorization of 
effective stiffness matrix is necessary for obtaining the solution. In this work, an implicit 
method called Newmark method is used for integration. 

The Newmark method makes use of the following assumptions [4]: 

t+At {£/} = *{#} + [(1 - 6 ) *{£/} + S' t+At {U}] At 

t+A *{17} = t {U}+ *{l7}At+[(i- «)*{£/} + a t+At {U}\(At) 2 (3.31) 

The parameters a and 8 can be determined so as to obtain desired accuracy and stability. 
The constant average acceleration method corresponding to a = j and 5 = \ is used in 
this work, which results in an implicit unconditionally stable integration scheme. The 
Newmark constants are defined below: 


1 


Oo = 

a (At ) 2 


1 

Oi = 

a (At) 


1 

02 = 

2a ~ 

a 3 = 

1 

i — 1 

<1 


04 = 8 At (3.32) 

Substituting for t+At {C7} from equation (3.30) and splitting t+At {[/} as: 

t+At {U} = t {U}+ t {AU} (3.33) 

the expression (3.30) becomes: 

‘[K d ],{AU}+ (oo'(M]'M+ '{/»= ‘ +i W 


(3.34) 



where the effective stiffness matrix t: \K d \ and the effective force vector are given 

by: 

^Kd] = o 0 *[M] + *[. K ] (3.35) 

t+M {Fd] = t+At {F} + *[M]( oo £ {C/} + Q\ *{11} + g 2 *{U}) (3.36) 

Since the dynamic equilibrium is satisfied at time t and assuming that change in mass 
matrix in an increment is not very significant, the following equation must hold, which 
can be obtained from the finite element discretization of the integral (2.41) written at 
time t: 

*[M\ *{U} + *{/} = *{F} (3.37) 

Using equation (3.31), this can be expressed as: 

a 0 t [M] t {U}+ *{/}= *{F d } (3.38) 

Decomposing the effective force vector as: 

t+At {F d }.= *{F d }+ t {AF d } (3.39) 

and using the expression (3.38) for t {F d }, the following equation in t{AU} is obtained 
from equation (3.34): 

V K d] t{AU} = t {AF d } (3.40) 

Using equations (3.36) and (3.38) we get: 

t {AF d } = t+At {F d } - ‘{F d } = t+At {F} - *{/} + 'ft} + 02 *{U}) (3.41) 

This algebraic system (equation (3.40)) is solved to obtain the incremental displacement 
vector t {AU}. 

Once the incremental displacement is obtained by solving equation (3.40), the velocity 
and acceleration at time t + At are obtained using the following expressions: 

t+at {fr} = 0<,(‘ +A, {C/}- d! ■{(/}- 02 '{&} (3-42) 

«•**{(/} = ‘{V}+ a a ‘{U}+ a, HA, {U} ( 3 . 43 ) 

3.4 Determination of Stress 

The evaluation of stress components (at the Gauss points of the elements) is done by the 
following stepwise procedure: 



(3.44) 


1. Calculate the relative deformation gradient £ +At [F]: 

t+Atpi £ j 9(i(A«{) 

* ij ~ ij + ~d^~ 

It should be noted that the position vector x corresponds to the equilibrium position, 
though t{Au} may be incremental or iterative displacement (See section 3.9.3). 
If the reference position vector corresponds to a non-equilibrium position, say the 
configuration at t^ , there exists a possibility of divergence of the numerical method. 

2. Decompose t+At [F] = t+At [J?] t+At [U] and determine t+A | [Q] and t+A t[U p ] using 
equation (2.3). 

3. Determine the principal incremental logaritmic strain t Ae L using equation (2.15). 

4. Calculate t {Acr M } using equation ( 2.39). The integration of the constitutive equa- 
tion is performed using the Euler forward technique described below. 

5. Use equation (2.12) to calculate the Cauchy stress components at time t + At 


3.4.1 Integration of the Constitutive Equation 

Various techniques exist for the integration of the constitutive equation. A simple and 
robust technique is the Euler forward integration scheme discussed below. 

Assume that the principal incremental strain components have been calculated and 
the state of the Gauss point at time t (elastic or plastic ) is known. 

• If the state at time t is elastic, 


1. Calculate the stress increment assuming elastic behaviour. 

t {Aa M } = (C E '] t {Ae L } 


(3.45) 


2. Calculate t+At {cr} using (2.12). 

3. Determine V eg and t+At a eq using equation (2.19). 

4. If t+At < 7 e , < t c y , then the Gauss point is still in elastic state, or if t_ ^ A V e g 
t a eg the Gauss point is neutrally loaded. Return. 

Else, a transition from elastic to plastic has ocurred. Calculate: 


t+At- 


Ratio = 


; eq 


t+A^j 


(3.46) 


eq 


'eq 



and change the state to plastic . The sub-incrementation method is followed 
and the stress components with respect to the material frame are updated after 
each subincrement by the increment in stress components corresponding to the 
elasto-plastic strain sub-increment. The [C EP '] corresponding to last updated 
state is used : 


,+A V*}® =‘ +a ‘ +•+“ dt{Ae L } fori = l,n (3.47) 


where, 

L = (1 - Ratio), {At 1 } 
n 

l+A!{ ff M}(0) = t + Ratm [ C E'J t { Ae I } 

t+4l[ C EP](°) _( [ C EP] eva l mital at e,Atf ff Mjrai 


(3.48) 

(3.49) 

(3.50) 


Use equation (2.12) to find t+At <7y. Return. 


• If the state at time t is plastic, check for unloading as described in section (3.4.2). 
If there is no unloading, the sub-incrementation method described in (4) above is 
applied with Ratio set to zero. 

This describes the Euler forward integration scheme for the elasto-plastic model de- 
scribed in section (2.7) 


3.4.2 Unloading Scheme 

The phenomenon of local unloading is to be incorporated to reproduce more closely the 
elasto-plastic response of a structure. The unloading criterion defined by Chakrabarty [14]: 

*{n} T t {Ae T }<0 (3.51) 

where, t {Ae} L is the incremental strain and *{n} T represents the unit outward normal 
to the yield surface at the current stress point {V} in a 9-D stress space, is imple- 
mented in this work. Neutral loading and postive loading are represented by replacing 
the “<” symbol in( 3.51) by “=” and “>” symbols respectively. 

When unloading is detected, the yield stress at that Gauss point is changed to the 
equivalent stress at the Gauss point at time t (f<j y = t o’ e q), the state tag of the Gauss 
point is changed from plastic to elastic and the elastic constitutive equation is applied to 
calculate the incremental stress, ‘[Acr^]. 



3.5 Skyline Scheme of Global Assembly 


The three-dimensional dynamic elasto-plastic contact problem is computationally inten- 
sive one and therefore every effort should be made to reduce the storage and the compu- 
tational time required for solving the problem. Skyline assembly scheme [4] is one such 
technique, which significantly reduces the storage required for the global stiffness ma- 
trix ‘[If]. In this storage scheme, only the elements below the skyline of the matrix ‘[if] 
are stored in a one-dimensional array t {A}. However, along with the storage scheme, a 
specific procedure for addressing the elements of ‘[If] in t {A} is needed. Thus, before 
proceeding with the assembly of the element stiffness matrices, it is necessary to establish 
the address of the global stiffness matrix elements in one-dimensional array ‘ {.A}. 

The element pattern of a typical global stiffness matrix is shown below: 


*[K} = 


ku k\2 0 ku 0 

k\2 k%2 k‘23 0 0 

0 k‘23 ^33 &34 0 

ku o k 3 4 k 44 k 45 

0 0 0 &45 &55 


Since the matrix ‘[If] is symmetric, only the elements above the diagonal and the diagonal 
elements need to be stored. Defining by m* the row number of the first non-zero element 
in column i, the variables mi = represent the skyline of the matrix, and the 

value (m i+ i — mi) represents the column height. (The zero elements below the skyline are 
stored, since they may become non-zero during the solution process). The column heights 
are determined from the connectivity arrays of the elements. Once the column heights 
of the global stiffness matrix are known, we can store all the elements below the skyline 
of ‘[If] as a one-dimensional array in ‘{A}; i.e., the active elements (elements below the 
skyline) of ‘[If] are stored consecutively in ‘{A}: 


^-fA}^ — {Ain, ^“22) ki2, k 33) A)23j k 44 ] &34j 0, ^14; k 33 , k 43 } 


In addition to ‘{A}, we also define an array {MAXA}, which stores the addresses of 
diagonal elements of ‘[If] in ‘{A}; i.e., the address of the i th diagonal element of *[K\, ku, 
in array ‘{A} is given by MAXAi. Thus, 


{. MAXA} t = { 1 , 2 , 4 , 6 , 10 , 12 } 



Note that, MAXAi is equal to the sum, of column heights upto the (i — 1 ) t/l column plus 1. 
Hence the number of active elements in the i th column of is equal to MAXA i+l - 
MAXAi , and the element addresses are MAXAi, {MAXAi + 1), . . (MAXA i+l - 1). 
It follows that using the above storage scheme of t [K\ in 4 {,4} . together with the address 
array {MAX A} as defined above, each element of 4 [jK] in ‘{A} can be addressed uniquely. 

3.6 Static Condensation Scheme 

The contact analysis is done in an iterative manner to determine the nodes which are 
out of contact, in sticking contact or in slipping contact. Solving the complete system of 
equations in each iteration is highly inefficient in terms of time and computer resources, 
especially when the contact region is small compared to the whole system. In order to 
reduce the solution time, it is desirable to condense the effective stiffness matrix and the 
effective force vector to the size of the degrees of freedom to which the contact boundary 
conditions are to be applied. 

The condensation procedure starts by dividing the nodes into two categories: the 
nodes at which the contact boundary conditions are to be applied (referred to as the 
nodes of type 1 ) and the remaining nodes (referred to as the nodes of type 2). To facilitate 
the condensation, the nodes are renumbered so that the nodes of type 1 are numbered 
first (called internal node numbers, while the original node numbers are called external 
node numbers ). Then the essential boundary conditions (i.e., the displacement boundary 
conditions) are applied to equation (3.40). 

Then these equations are partitioned as follows: 1 

’ [K n ] [Kn] 

_ [K 2l ] [k 22 ] 

where {Ai7i} denotes the incremental displacement vector corresponding to the nodes 
of type 1 and {A U 2 } denotes the incremental displacement vector corresponding to the 
nodes of type 2. Equation (3.52) can be seperated as follows: 

[K n ] {Ai7i} + [K n ] {AU 2 } = {AFJ (3.53) 

[K 2l ]{AU 1 }+ [K 22 ]{AU 2 } = {A F 2 } (3.54) 


< 

{AC/,} 

y zzz < 

{AF,} ^ 


l {Aty J 


[ {AFi} 


(3.52) 


x In this section, the superscript/subscript denoting time have been omitted for the sake of convenience 



Solving for {A U 2 } from equation (3.54), we get: 


{AU 2 } = [K 22 ]~ l ({A F 2 } - [K 2l ] {AU,}) (3.55) 

Substitution of this expression for {AC/ 2 } in equation (3.53) and rearrangement of the 
resulting equation leads to the following condensed set of dynamic equilibrium equations: 

[K]{AU} = {A F} (3.56) 

where, 

[K] = [K n ]-[K 12 ][A] (3.57) 

{AC/} = {AUi} (3.58) 

{A F} = {AFi} — [K 12 ] {B} .(3.59) 

[A] = [K 22 )~ l [K 2l ] (3.60) 

{B} = [K 22 ]~ l {AF 2 } (3.61) 

The number of equations in the condensed set (3.56) is equal to the number of degrees 
of freedom of type 1. The coefficient matrix [K] and the right side vector {A F] are 

evaluated from equations (3.57 - 3.61) using the partitioned matrices and vectors of 

equation (3.52). Note that while evaluating [A] and {B}, it is not necessary to invert the 
matrix [K 22 ]~ l . Instead one can solve the following equations by the Gauss elimination 
method: 

[K 22 \ [A] = [K n ] (3.62) 

[K 22 ]{B} = {AF 2 } (3.63) 

To reduce the time further, one can store the upper triangular form of [K 22 ]. In 
this way, one can solve for the columns of [A] by performing Gauss elimination and 
back substitution operations on the corresponding columns of [/Gi]- The vector {//} is 
obtained in similar fashion by performing the Gauss elimination and back substitution on 
the vector {AF 2 }. 

3.7 Modified Newton-Raphson Scheme 

The solution to equation (3.40) represents only an approximate solution to the governing 
equation (2.44), because of the linearization and simplifications employed in deriving the 



equation (2.o2). The error obtained should to be minmized and there are various iterative 
techniques for this. One of the simplest and very effective technique is the Modified 
Newton- Raphson algorithm, which offers fast convergence with less computation [4]. 

This technique can elegantly be represented as follows: 

Solve: 

t{AU}M= t {AF u }^ (3.64) 

where, 

,{A F„}V-» = (3 65) 

= ,{AF ( } (3.66) 


till the convergence criterion: 




< toL 


||t + At{ i r}(i)|| 

is satisfied. The vector t{AjF u } is called the unbalance force vector. 


(3.67) 


3.8 Divergence Handling Procedures 

The modified Newton Raphson method diverges in some cases; i.e., the ratio of unbalance 
force to the external force, instead of reducing, increases in each iteration. Also, in 
some other cases, the rate of convergence is not fast enough, requiring a large number 
of iterations for reducing the unbalance force to acceptable limits. Two of the methods, 
which are simple and fairly effective for handling divergence or accelerating the rate of 
convergence are Under-relaxation method and Line search method [8] . These methods 
are generally applied by modifying the displacement vector t {Au}^. However, in the 
contact analysis scenario, it is not acceptable, since the contact status may get changed 
due to modification in the displacement vector. Therefore, in this work, these methods are 
applied by modifying the unbalance force vector ^Ai^p -1 ) and repeating the solution 
process when divergence is detected. 

In this work, the Newton Raphson iterations are said to be diverging when: 

. J,{AF„}<‘- 1 >|| 

ni <+a w <w) ii 

is satisfied. The value of the factor fi is taken as 0.9. 

The under-relaxation and line search methods are discussed below: 



ll.{AK} (i) j| 

||<+At{ f }(i)|| 



Under-relaxation method: If divergence is detected in the i th iteration, the unbalance 
force vector t{AF u }^ ^ is scaled by a factor a repeatedly, till the above criterion is 
not satisfied: 

t{AF u }^= a t {A F U }W (3.69) 

If the under-relaxation is not able to find a value of unbalance force for which 
the above criterion is not satisfied, the solution run is terminated with a warning 
message. 

Line search method: In this method, the basic aim is to search for the point, which 
corresponds to minimum unbalance, along a direction specified by the unbalance 
force vector t {AF u }^ l \ However, it is not very efficient to perform a full line search, 
since, in that case the number of solution runs required may become prohibitively 
high. The approach used is to repeat the solutions for a fixed number of values of 
factor aj (defined below), within a preselected range. 

t{AF«} (i_1) = a# t{AF u } (i-1) (3.70) 

where, a min < a j < a mai 

The value of aj corresponding to the minimum unbalance force is chosen as the 
best point and the solution corresponding to this o j is accepted as the result of line 
search method. 

Since in both the above methods, the complete iteration is repeated, including a 
number of contact iterations, it is very costly to do either under-relaxation or line search, 
in terms of both computational time and resources. Therefore, it is preferable to use 
these methods only when some difficulties in convergence are observed. Under-relaxation 
is cheaper compared to line search, while line search is more efficient. 

3.9 Numerical Aspects 

The selection of parameters controlling the numerical scheme in dynamic non-linear analy- 
sis is of importance not only for efficiently obtaining the solution, but also for the accuracy 
of the solution. The following sections discuss various aspects related to the numerical 


scheme used. 



3.9.1 Choice of Time-step 


A procedure for the selection of time step is given in Bathe [4], which may be used to select 
an initial time step and mesh. A time, step and mesh convergence study should then be 
performed to arrive at the proper values. It is important to note that an unconditionally 
stable method does not provide a solution close to the exact solution for an arbitrary time 
step. Phenomena such as amplitude reduction and period elongation may occur due to 
an improper choice of time step. Hence the proper way to perform a dynamic non- lin ear 
analysis is to study the effect of time step and mesh size on the solution and to search 
appropriate values. 

The choice of a time step is based on the type of problem considered. In general, 
dynamic problems can be categorized as: 1. Strucural dynamics problems and 2. Wave 
propagation problems. 

Structural Dynamics Problem: The step by step procedure for modeling a strucutu- 
ral dynamics problem is given below. 

1. Identify the frequencies involved in loading, using a Fourier analysis if nec- 
essary. These frequencies may change with time. Let the highest significant 
frequency contained in the loading be w u . 

2. Choose a finite element mesh that can accurately represent the static response 
and accurately represent all frequencies upto about w co = 4w u , where Wco is 
called as the optimum critical frequency. 

3. Perform the direct integration analysis. The time step At for this solution 
should be equal to about where T co = 

Wave Propagation Problems: Assuming that the critical wavelength to be calculated 
correctly is L w , then the time taken by the wave to travel past a point is given by: 

u = — (3.71) 

c 

where c is the wave speed. Assuming that n time steps are necessary to represent 
the travel of the wave, the time step should be: 

n 


(3.72) 



and the effective length 2 of a finite element should be: 

L e = cAt (3.73) 

The above choice of time step and corresponding effective length represents the 
complete wave travel accurately. But, care has to be taken to ens ure that At < — , 
where T n is the smallest time period of the problem. 

It is found that the most accurate solution is obtained by integrating with a time 
step equal to the above limits (denoted by At cr ) and the solution is less accurate, when 
a smaller time step is employed. This deterioration in the accuracy of the predicted 
solution when At is smaller than A t^ is more pronounced, when a relatively coarse 
spatial discretization is used. 

3.9.2 Deformation Dependent Loading 

Deformation dependent loading is taken care of by evaluating the consistent force vec- 
tor t+At {F} at time t + At for the calculation of the unbalance force vector (3.65). An 
alternate procedure given by Bathe et al. [17] incorporates deformation dependency into 
the formulation of the governing equation, but this results in an asymmetric stiffness ma- 
trix, the solution of which is computationally expensive. The method used in this work 
may take a slightly more number of iterations for convergence, but this loss is offset by 
the gains acquired by retaining a symmetic stiffness matrix. 

3.9.3 Stress Updation 

There exist two methods for updating the stress at a Gauss point: 1. Iterative wpdation 
and 2. Incremental updation (Crisfield [8]). 

Iterative Updation: The iterative displacement is used to compute the iterative strain, 
which is then used in the constitutive equation to update the stress at iteration (i— 1) 
to stress at iteration (i). 

Incremental Updation: All the iterative displacements upto iteration (i) are added to 
find the cumulative incremental displacement of the current increment. This is used 
to compute the incremental strain, which is used in the constitutive equation to 


2 Effective length is the smallest distance between any two nodes in the mesh used 



update the stress from last equilibrium point (end of the previous increment) to the 
state after iteration (*). 

The first method may lead to spurious unloading, since the unbalance force vector 
calculated in some iterations may be in different directions compared to that of the initial 
force vector. This can be avoided by using the incremental updation procedure. The 
examples solved in this work employ the incremental updation method. Iterative updation 
may be used if the convergence is monotonic. 



Chapter 4 


Contact Formulation 

The dynamic, large deformation, elasto-plastic, updated Lagrangian finite el em ent for- 
mulation developed in the previous two chapters, can be used for contact analysis with 
appropriate modifications. The most important modification is the development of an 
additional set of finite element equations (involving the contact stiffness matrix) relating 
the unknown contact forces and displacements. These equations, which consist of the 
kinematic constraints and the contact force expressions based on the discretization of the 
contact boundary, are developed using the Lagrange multiplier method. These develop- 
ments are described in the following sections. Finally the algorithm for the analysis of 
dynamic, large deformation, elasto-plastic contact problem including the effects of friction, 
is described. 


4.1 Contact Constraints 

A typical two body contact system is shown in Figure (4.1). The bodies occupy the 
domains t+A A fl and t+A V 2 at time t + At. The boundaries of t+A ^ 1 and t+Af V 2 are 
denoted by t+At S l and t+At S 2 and their interior volumes by t+A V 1 and t+A V 2 respectively. 
At any time t + At, the boundary of contact body n (n = 1,2) can be partitioned as: 

t+Mgn _ t+Atgn (J t+Atgn y t+Atgn (41) 


where, 

t+At s £ = surface on which displacements are specified, 
t+At 5 n _ sur f ace 0 n which tractions axe specified, 
t+At S n _ surface on which contact occurs. 




Figure 4.1: Two body contact system 

Assuming the boundary t+At S n to be smooth, outward unit normal vector at point 1 t+At x n 
on the boundary t+At S n is denoted by t+A W". Other two orthonormal boundary vectors 
(in tangential directions) are t+At N£ and t+A W£ (see Figure (4.2)). 

Suppose that two boundary points t+At x l and t+At x 2 are in contact with each other at 
time t + At and the contact traction at t+A V* is denoted as t+At T£. Then by Newton’s 


third law of motion, we have: 


t+Atrjil = ___ t+Atrp2 


where t+At T£ can be expressed as: 

3 

t+Atrpn __ ^2 t+At^n t+Atpjn ^ 

t=l 

where, t+At i- 1 denotes the component of t + At T^ along *+ At jV? . 

If the two contact bodies are not. welded together, then the normal component of 
traction can not be tensile. Thus we have: 


t+At+n 


t? < o 


Further, from the Coulomb friction law we have: 


\t+At t ny + (t+At^)2< 


1 Vectors are denoted by a letter with underbar 



Hitting Body (Body 1) 



Figure 4.2: Points in contact and associated normal and tangent vectors 

In equation (4.5), the “<” symbol represents the sticking friction condition and ’ 
represents the slipping condition. Further, the resultant tangential (friction) force is in 
the direction opposite to the relative velocity of slipping. The equations (4.4) and (4.5) 
represent traction constraints. 

Physical considerations require that one body can not penetrate the other. Thus, for 
the points in contact, pentration must be zero. This condition can be mathematically 
stated as: 

*+ A ^ x = ( t+M x 2 - t+At x l ). t+At Nl = 0 (4.6) 

where t+At p x represents the penetration between the two points in contact. In addition 
to (4.6), we need to have two more conditions to enforce the sticking friction condition, 
which ensure that there are no relative displacements (slip) between the contacting points 
in the tangential directions. These conditions can be stated as: 

t + A * p 2 = (*+ A ^ 2 - t + z V ). t+At Nl = 0 

t+A1p 3 _ (t+A^2 _ t+At^ly _ Q 


The conditions (4.6), (4.7) and (4.8) are called as the kinematic constraints. 


(4.7) 

(4.8) 



4.2 Contact Force Expressions 


As the contacting bodies are discretized for the purpose of developing finite element 
equations, the contact boundary gets automatically divided into surface elements, called 
as contact segments. The contact force formulation in this work employs a node-to-surface 
interface model, since this work deals with a large deformation problem including slipping 
at the contact interface. Normally in a contact formulation, contact is considered only at 
discrete nodes. Here, it is assumed that a node of the hitting body (body 1) comes into 
contact with a segment of the target body (body 2). The first step in the development of 
the contact force expression is the determination of the unit normal and tangent vectors 
at a point on the target segment. This is described in the next paragraph. 

Geometry of a target segment can be described using 2-D shape functions. Thus, 

‘ +i k 2 (f,o) = E%K^) ,+i k? («) 

j= 1 

where, <&_,(£, g) j = l,n are the shape functions, t+A< x| denotes the position vector of 
node j and N is the number of nodes per target segment. The side of a segment on which 
contact can occur is referred to as the positive side of the segment and the other side as 
the negative side of the segment. The local numbering of the nodes on a contact segment 
should be counter clockwise when we face the positive side of the segment. This is to 
get the correct sense for the outward normal vector at a given point. A normal vector at 
point t+A *r 2 (£, v) on the positive side of the segment is defined as: 


t+Atjy2 _ 

. t+Atjy2 ^ £+Atjy2 

(4.10) 

where, t+At iV| and t+A W 2 are the two tangent vectors at t+A ^ 2 (£,7?) given by: 


t+Atjy2 

= ‘ +i k 2 K,»), f 

(4.11) 

t+At iv; 

fr 

's? 

f 

t 

II 

(4.12) 

The unit normal and tangent vectors are defined as: 



_ t+Atjy2/|t+Atjy2| 

(4.13) 

t+At N_l 

= t+At N^/\ t+At N^\ 

(4.14) 

HAt N% 

_ t+ Atjy2 j j t+At^j-2 j 

(4.18) 


Since the contact is considered only at discrete hitting nodes, we cam deal only with 
point forces at hitting nodes rather than with contact tractions. At time t + At, let 



t+At{ F 2 } be an array of Cartesian components of point force vector t+At F 2 at a point of 
the target segment at which the hitting node makes contact. By Newton’s third law, the 
force vector on the hitting node will be — t+At {F 2 }. 

Let t+At {<5u 2 } and t+At [d i/ 1 } be the virtual displacement vectors at time t + At at 
these points (on the target and hitting bodies respectively), when they come in contact. 
Then, the virtual work at hitting node can be expressed as: 2 


5w c = t+At {5u 2 ~ SuY^F 2 } (4.16) 

Note that the elements of t+At {(5« 1 }, namely 5 t+At u 1 , S t+At v l , 8 t+At w l are the virtual 
nodal displacements of the hitting node. The elements of t+At {<J m 2 } can be expressed in 
terms of the virtual nodal displacements of the target segment (8 t+Ai uf, 8 t+At vf, 5 t+At wj, 
i = 1, N) using the 2-D shape unctions ($j, i = 1, N). Thus: 


t+Af {^ 2 - 5u 1 }= [Q c ] t+At {«5 Wc } 


(4.17) 


where, 


[Qc] — 


-10 0 0 0 $2 o 0 ... $ N 0 0 

0-10 0 0 0 $ 2 o . . . 0 o 

0 0 —1 0 0 $1 0 0 $2 • • • 0 0 


(4.18) 


and 


t+At {* u c } T = {5 t+A V, t+A V, < 5 t+Ai w l ,5 t+A W v 8 t+A V u 5 t+A hl . . . , S t+A H> (4-19) 


4.2.1 Sticking Friction Condition 


For the sticking friction, the vector t+At {F 2 } can be expressed as: 



t + A W 1 2 1 t+At N%x t+At lV | 1 


f ^ 

t+Atf2 


t+At N? 2 t+At N%2 t+At /Vf 2 

< 

t-\-At^ 2 


t + A W 2 3 t+At N 2 z t+At /V | 3 


t+Atf2 


(4.20) 


or, 


t+At{ F 2^ _ (4.21) 

where, t+At /Vj represents the j th Cartesian component of the unit vector t+A W 2 and t+At f? 
is the component, with respect to t+A 7V 2 , of the contact force vector t+ t F. - Combining 


2 In this section, the right superscript denoting the iteration number has been omitted for convenience 



equations (4.16), (4.17) and (4.21), we get the following expression for the virtual work 
at the hitting node: 


6w c = HAt {5u c } T t+At {r c } 

(4.22) 

where the vector t+At {r c } is given by: 


t+At {r c } = [Q c ] T 4 +At [lV 2 ] t+At {f 2 } 

(4.23) 

This is the contact force expression at a sticking node (hitting node, 
sticking friction condition). 

which is under 

4.2.2 Slipping Condition 


The contact force expression for a slipping node (hitting node, which is slipping on the 

target segment) is developed below: 

The vector t+At {r c } in equation (4.23) can be split as a contribution from the normal 
contact forces and tangential contact forces. Thus, 

1+A ‘{rJ = ,+A ‘K} + ,+A ‘{r.} 

(4.24) 

where, 

t+At {r n } = [Q c ] t t+At {N 2 } t+At f 2 

(4.25) 

and 

t+At {r t } = [Qcf ( t+At {lV 2 2 } i+Ai /| + t+At {N 2 } t+At fi) 

(4.26) 

Here, the arrays t+At {N?}, i = 1,2,3 contain the Cartesian components of the vec- 


tors t+A t /Vf . 

The three components of a contact force are not independent at a slipping node. Prom 
Coulomb’s friction law, we have the relation: 

v /(t+At / 2)2 + ( t+ A t / 2)2 = ^|*+A‘/2| (4.27) 

where n is the coefficient of friction at a slipping node. Therefore, the tangential compo- 
nents of the contact force can be written as: 

t+At /f = ~ t+At /i ^ cos 0 

t+Atf2 = _t+Atj2^ sind 


(4.28) 

(4.29) 



where, 6 is the angle between the resultant tangential (friction) force and the vector 
and since the normal contact force is always compressive, ” sign is used in the right 
hand side of both the above equations. 

Calculation of the angle 6: 

The friction force in the case of slipping is in a direction opposite to the relative 
velocity of slipping, which is same as the difference in incremental displacement of the 
target point and the hitting node. Therefore, the direction of the friction force on the 
target body is given by: — {*Au 2 - tAu 1 }* = — ([Q c ] t{Au c }) t , where the right subscript t 
indicates the tangential component. Define: 

,Ah 2 = - ,+i! {«f} r [Q c ] ,{Av c ] (4.30) 

t A«, = - ,{&»„} (4.31) 

Thus, the angle 9 is given by: 

“‘(ill) < 4 - 32 > 

However, to avoid the difficulty in handling infinity, the computer implementation uses 
cosine function instead of tangent function for calculating 6. 

Substituting equations (4.28) and (4.29) in equation (4.26) and then combining equa- 
tions (4.25) and (4.26) with equation (4.24) we get: 

t+At {r c } = [Q c f ( t+At {iV?} - t+At {iVf} n cos 9 — t+At {iVf } /x sin O) t+At fi (4.33) 

4.3 Kinematic Constraints in Nodal Form 

The kinematic constraints should be written separately for the sticking nodes and for 
slipping nodes. 

4.3.1 Sticking Friction Condition 

To express the kinematic constraints in terms of the nodal displacements, equations (4.6), 
(4.7) and (4.8) are combined in matrix form for the i th Newton-Raphson iteration as: 

t+At|p|(i) _ t+ t+At {x 2 - x 1 }® (4.34) 

where, 


= {'+^ 5 °. 


(4.35) 



We have, 

t+At {x 2 - = t+A \x 2 - + t {A u 2 - A u 1 }® (4.36) 

where, the second term on the right hand side is an array of Cartesian components of 
the difference of the displacement vectors of the target and hitting bodies at the contact 
node, for the current iteration. For the first iteration: 


t+At {x 2 - x 1 }^ = « At {x 2 - = V - x 1 } 


(4.37) 


Analogus to equation (4.17), the second term of the right side of equation (4.36) can be 
written as: 

t{Au 2 - Au 1 }® = [Q c ] *{Au c }« (4.38) 

where, the vector: 


{Au c } (l)7 = {*Au lW , t Au lW , t Au; lW , t A uf l \ tAu 2W , t Au> 2(l) , . . . , t Aio^ W } (4.39) 


contains the nodal values of displacement at the hitting node and the nodes of the target 
segment for the current iteration. Recognizing the first term of the right hand side of 
equation (4.36) as corresponding to (penetration for the previous iteration), 

equation (4.34) can be written as: 

t+At {p}« = t+At {p yi-i) + t+At^Y [Qc] (4.40) 

where for the first iteration, 

= t+At^(o) = (4. 41 ) 


Since t+At {p} (i) is zero, the kinematic constraint at a sticking node is: 

t + At {p y(i-i) + t + At^ T [Q c ] t {Au c }« = {0} (4.42) 

The first term of the right hand side of equation (4.36), on similar development, gives: 

t+A^i-r) = t+At^Y (4.43) 

where, t+At {x c } (i-1) is the array of coordinates of the hitting node and the nodes of the 
target segment updated upto the previous iteration. Expression (4.43) is used to calculate 

t+At {p} {i ~ 



4.3.2 Slipping Condition 


Development of the nodal constraints for the slipping nodes is similar to the sticking 
nodes, except that we have only one kinematic constraint, viz., equation (4.6). 

We can write the matrix form of the kinematic constraints for the i th Newton-Raphson 
iteration as: 

t+At p? = t+At {A r 1 2 } T t+M {x 2 - x 1 }® (4.44) 

Substituting equations (4.36) and (4.38) in equation (4.44) we get: 

+ ‘ +i, {ivy [QJ (4.45) 

where for the first iteration, 

t+Atpii-D = t+ a^(0) = tpi ( 4 .46) 

Since t+At p^ is zero, the kinematic constraint at a slipping node is: 

t+A^i-i) + hm {N 2 } t [Qc] t{A Ue }(i) = o (4.47) 

Similar to the sticking condition, we have: 

t+ A ^i-D = t+At^r jqj (4 48) 

where, t+At {x c }^ -1 ^ contains the coordinates of the hitting node and the nodes of tar- 
get segment updated upto the previous iteration. Expression (4.48) is used to calculate 


4.4 Lagrange Multiplier Method 

In Lagrange Multiplier Method, the contact forces are considered as primary unknowns 
and the kinematic constraints are enforced exactly. We have the contact force expression 
for a sticking node given by equation (4.23) and for a slipping node by equation (4.33). 
The kinematic constraint for a sticking node is given by equation (4.42): 

t+At [iV 2 ] T [Q c ] t {Au c } (i) = - t+At {p} (i_1) ( 449 ) 

and for a slipping node is given by equation (4.47): 

t+At{ N if [Qc] t { Au c } (i) = - t+Ai pt l] 


(4.50) 



Combining equations (4.23) and (4.49) for a sticking node, we get 


0 - " 

i 

l ‘+ A ‘{f}i0 } 

f t+At {p} (£_1) 

* +A, [«i] 0 


{ t{A“J (i) J 

\ - *+ A *{ rc }(i-D 


where, 

,+A %if = ,+a, [iV 2 ] T [Q c ] 

Similarly, combining the equations (4.33 and 4.50) for a slipping node, we get 

*+V i_1) 1 

- t+At {r c } (i_1) J 

where, 


0 

t+At te} 


t+At {<hY 

0 




(4.51) 


(4.52) 


(4.53) 


,+A ‘{«3 } T = ‘ +A1 K} r »J («4) 

•+ a, { ?2 } = [Q,f (‘-''“{iV, 2 } - ,+A, {f^} /i cos 0 — 1+i ‘{jV(} fi sin D) (4.55) 


Assembling these equations (4.51) and (4.53) over all the potential contact nodes, we 
get the following global equation: 

0 - »*[Q x f 

_ - t+At [Q 2 ] 0 

Here, the vectors t {A/7 c }^^ and are the assembled versions of t{Au c }W and 

£+At {r c }0— 1 ) respectively. Further, the vectors t+At {F 2 }^^ and t+At {P}( i_1 ) are the assem- 
bled versions 0 f t+At {/ 2 }(* ) and t+Ai {p} (i-1) respectively for the sticking nodes and t+At / x 2W 
and for the slipping nodes. Similarly, the matrices t+A \Qi} T and t+A1 [(3 2 ] axe 

global versions of t+At {<7i} r and t+At {gi} respectively for the sticking nodes and t+Ai {g3} T 
and £+At {g 2 } for the slipping nodes. 

Note that, while the whole global coefficient matrix is known from the geometry, only a 
part of the right side vector, namely t+At {P is known from the geometry. The other 
part, t+At {P c } (i_1) , is unknown. This vector is eliminated by combining this equation 
with equation (3.64) of chapter 3. Before this, the set of equation (3.64) is condensed 
to retain only those degrees of freedoms, which axe associated with the contact nodes. 
This condensed and combined set is solved to find the nodal displacements of the contact 
nodes and the contact forces. Then equation (3.55) is used to find displacements of the 
nodes of type 2. 


t+Atj- p2^(i) 

t{AE/c}« 


= < 


t+At^pyi-l) 

- t+At {R c }^- l ) 


(4.56) 



4.5 Algorithm used for Dynamic Large Deformation 
Elasto-plastic Contact Problem 

For solving a typical dynamic, large deformation, elasto-plastic, contact problem, the 
following steps are used: 

1- Finding the potential contact nodes for which the hitting and target nodes are closer 
than a prescribed length. 

2. Renumbering the nodes: Nodes of the hitting and target bodies are renumbered 
such that the contact nodes are numbered first. This is carried out to facilitate the 
static condensation of the stiffness matrix which helps in reducing the computational 
time. 


3. Forming the coefficient matrix and right side vector: Based on the current geometry 
and state of stress, condensed form of effective stiffness matrix and force vector are 
formed. 


4. Beginning the Newton-Raphson iteration. 

5. Contact search: Search for the target segment corresponding to each hitting node. 
Master-slave algorithm is used here. 

6. Contact iterations: 


• Initially all the potential nodes are assumed to be in sticking friction condition. 
Then, the contact stiffness matrix and right side vector are formed and com- 
bined with the condensed form of the effective stiffness matrix and the effective 
force vector. This system of equations is solved. 

• Then, the “out of contact” nodes, for which the normal component of the 
contact reaction is tensile (refer to equation (4.4)) are found. The contact 
reactions at these nodes are set to zero during the subsequent contact iterations. 


• After removing all the “out of contact” nodes, the nodes which are slipping are 
determined. We know that a node slips, if the contact reactions at that node 
violate equation (4.5), expressed in a nodal form. 


yjsftvP? W"." 
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• The contact iterations are repeated till the correct direction of friction force is 

obtained for all the slipping nodes and the status of all the contact nodes does 
not change. 

7. Finding the displacements of the non-contact nodes (type 2 nodes) from the dis- 
placements of type 1 nodes. 

8. Euler forward integration scheme: Depending on the choice of Stress measure (i.e 
Jaumann or New objective) and type of updation (i.e iterative or incremental), the 
Euler forward integration is performed. 

• J aumann Stress Measure: From the iterative or incremental displacement, first 
the strain is found. Then, the change of state (i.e elastic to elastic , elastic to 
plastic, plastic to plastic or plastic to elastic) is determined. Finally, depending 
on the change of state, Euler forward integration is performed. 

• New Objective Stress Measure: Using the iterative or incremental displace- 
ment, first deformation gradient matrix is found out and then using the polar 
decomposition theorem, the logarithmic strain is found. The initial stress is 
transformed to the material coordinate system. The change of state at the 
Gauss point is found. Depending on the change of state at the Gauss point, 
the Euler forward integration scheme is implemented. Finally, the updated 
stress is transformed back to the fixed frame. 

9. Updating: The stresses, strains and contact forces are updated. 

10. Convergence: The unbalance force is found and convergence is checked. If the 
Newton-Raphson iterations converge, the results at the end of the increment are 
stored and the next time increment is started. If the convergence criterion is not 
satisfied, the Newton-Raphson iterations are continued till convergence. 



Chapter 5 


Results and Discussion 

A finite element code for the analysis of 3-D, dynamic, large deformation, elasto-plastic 
contact problem with friction is developed based on the formulations developed in the 
chapters 2, 3 and 4. The stiffness matrix, mass matrix and force vectors are evaluated 
using numerical integration, while Newmark’s scheme is used to update the velocity and 
acceleration after each time-step. The code is developed for 8-noded brick element and 
2-point Gauss-Legendre numerical integration is used in each direction. 

The code with Jaumann stress rate and incremental Green-Lagrange strain measures 
was developed and validated by Jayadeep [57] for large deformation, dyn ami c, elasto- 
plastic analysis (for both the single-body as well as contact problems). Now, an additional 
set of objective stress and non-linear strain measures, namely new incremental objective 
stress tensor and incremental logarithmic strain tensor, is added. The objective of the 
thesis is to carry out the comparitive study of these two sets of stress and strain measures. 
This is done for the following four cases: 

1. Large Deformation, Dynamic, Elastic, Single Body Analysis. 

2. Large Deformation, Elastic, Contact-Impact Analysis. 

3. Large Deformation, Dynamic, Elasto-plastic, Single Body Analysis. 

4. Large Deformation, Elasto-plastic, Contact Analysis. 



5.1 Large Deformation Dynamic Elastic Single Body 
Analysis 

5.1.1 Dynamic Response of a Cantilever 


Plb/in2 

2 .& 



Time(sec) 


T, 


.X 


Figure 5.1: Cantilever : geometry and boundary conditions 


A benchmark problem solved by Bathe et al. [IT] has been chosen for the large de- 
formation dynamic elastic single body analysis. Figure (5.1) shows the geometry and 
boundary conditions of the cantilever analyzed. The geometric and material parameters 
used are: 

• Length of cantilever L = 10 inches 

• Height h = 1.0 inch 

• Width b = 1.0 inch 

• Young’s modulus E = 1.2 x 10 4 lb/in 2 

• Poisson’s ratio v = 0.2 

• Density p = 1.0 x 10 -4 lb-sec 2 /in 4 

Mesh and Time — step Convergence Studies 




A mesh and time-step convergence study is essential to ascertain that a proper mesh 
and time step are selected so that a correct solution is obtained with proper utilization 
of computational resources. First a mesh convergence study is carried out to select an 
appropriate mesh. As the problem is a dynamic one, so the time step convergence study 
is carried out next to select a proper time step. 

1. Mesh. Convergence: Following are the discretizations used for the mesh conver- 
gence study. The time step used is 1.35 x 10 -04 sec and the convergence tolerance 
of 1% is used. 

• Discretization No. 1: 

- Number of elements: 80 (20 x 2 x 2) 

- Degrees of freedom: 567 

• Discretization No. 2: 

- Number of elements: 640 (40 x 4 x 4) 

- Degrees of freedom: 3075 

From the mesh convergence study, it is found that the second discretization gives 
more accurate result. So, for the final study it is selected. 

2. Time — step Convergence: Following are the time-steps considered for the time- 
step convergence study. The discretization used is (40 x 4 x 4) and the convergence 
tolerance used is 1%. 

• Time-step 1= 1.35 x 10~ 04 sec. 

• Time-step 2= 0.45 x 10~ 04 sec. 

It is found that both the time steps give approximately the same results. For the 
final study, the time-step of 0.45 x 10 -04 is selected. 

The comparison between the Jaumann stress rate and new incremental objective stress 
measures is shown in Figures (5. 2-5. 3). The figures show that both the stress measures 
predict the same level of maximum displacement and equivalent stress for the single-body 
elastic analysis. 





5.2 Large Deformation Elastic Contact-Impact Anal- 
ysis. 

5.2.1 Impact of Two Identical Bars 
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Figure 5.4: Impact of two identical bars 

This problem is solved by Belytschko and Neal [44]. Figure (5.4) shows the initial 
configuration of the bars along with the boundary conditions. One of the bars is given 
an initial velocity so that it impacts with the second bar. The material properties are so 
chosen that the wave speed in each of the bar is 10.0 m/sec. Since c = the material 
properties chosen are: 

• Young’s modulus E = 10 N/m 2 

• Poission’s ratio v = 0.0 

• Density p = 0.1 Kg/m 3 . 

The geometric properties of each bar are as follows: 

• Length of the bar L = 10.0 m 


• Thickness = 1.0 m 



• Width = 1.0 m 


The initial velocity of the impacting bar i.e bar A is 1.0 m/sec. As Newmark’s scheme 
is used to calculate the velocities and accelerations, it has to be made sure that the hitting 
body travels with zero acceleration before the impact. In order to achieve this the initial 
gap between the two bars has been kept at 0.1 m and the contact search length has been 
chosen as 0.11 m. 

Mesh and Time — step Convergence Studies : 

First a mesh convergence study is carried out to select an appropriate mesh. Next, 
time-step convergence study is carried out to select a proper time step. 

1. Mesh Convergence: Following are the discretizations used for the mesh conver- 
gence study. The time step used is 0.1 sec and the convergence tolerance of 1% is 
used. 

• Discretization No. 1: 

- Number of elements: 10 (10 x 1 x 1) 

- Degrees of freedom: 132 

• Discretization No. 2: 

— Number of elements: 20 (20 x 1 x 1) 

- Degrees of freedom: 252 

The study shows that the second discretization gives better results. So it is chosen 
for further study. 

2. Time — step Convergence: Following are the time-steps considered for the time- 
step convergence study. The discretization used is (20 x 1 x 1) and the convergence 
tolerance used is 1%. 

• Time-step 1= 0.2 sec. 

• Time-step 2= 0.1 sec. 

The time-step convergence study shows that time-step 2 gives better results. So it 
is chosen for the further study. 



Figures (5.5) and (5.6) give the velocity time histories at the interface of first bar 
(Z = 10.1), as well as second bar (Z = 10.0). Figures (5.7) and (5.8) give the comparison 
between the maximum equivalent stress in bodies A and B respectively. Both the stress 
measures predict almost the same values of maximum equivalent stress. 












5.2.2 


Oblique Impact of Two Infinite Blocks 




x 


Figure 5.9: Oblique impact of two elastic blocks 

The secondproblem used for comparison is the oblique impact of two infinite elastic 
blocks (plane strain problem) shown in Figure 5.9. This problem is reported in refer- 
ence [49]. The top block is given a uniform initial velocity, Vo = (-10,0,-10) units (in 
the coordinate system shown in figure). The geometric and material properties are: 

• Hitting Block 

— length = 5 units 

- height = 4 units 

• Target Block 

- length = 16 units 

— height = 6 units 

• Young’s modulus ( E ) = 1000 units, for both the blocks. 

• Poisson’s ratio ( v ) = 0.0, for both the blocks. 

• Density (p) = 0.1 units, for both the blocks. 




Armero and Petocz [49] have solved this problem as a plane strain problem. Bhat [63] 
has solved the same problem, only for the frictionless case, using Jaumann stress rate 
measure and incremental Green-Lagrange strain measure, while Jayadeep [57] has solved 
the problem for both the friction and frictionless cases using Jaumann stress rate measure 
and incremental Green-Lagrange strain measure. In this work, however, new incremental 
objective stress and logarithmic strain measures are used and the results are compared 
to see the effect. Since the Poisson’s ratio is zero, the plane strain and plane stress 
formulations will give the same results. In this work, the problem is solved with only 
one element in the Indirection. To overcome the numerical difficulty arising due to low 
stiffness in F-direction, the displacement along this direction is fixed at all nodes for both 
the bodies. 

The following are the discretization and timestep considered: 

• Discretization: 

— Number of elements in the hitting block = 20 (5 x 1 x 4) 

— Number of elements in the target block = 96 (16 x 1 x 6) 

— Number of d.o.f in the hitting block = 180 
— Number of d.o.f in the target block = 714 

• Time-step = 0.01 units 

There are two different cases. In the first case, there is no friction at contact interface and 
in the second case, friction coefficient (/j) between the hitting and target bodies is 0.4. 

The X and Z displacements at point P are plotted against the time for both the 
Jaumann stress rate and new incremental objective stress measures in Figures (5.10) and 
(5.11). Both the displacements and time are measured from the instant of contact between 
the two blocks. Comparison of the maximum equivalent stress in bodies A and B is shown 
in Figures (5.12) and (5.13). The maximum equivalent stress predicted by the two stress 
measures differ at certain time-steps. No specific trend is observed. 
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igure 5.12: 
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Figure 5.13: Maximum equivalent stress in block B 





5.3 Large Deformation Dynamic Elasto-Plastic Sin- 
gle Body Analysis 

5.3.1 Dynamic Elasto-Plastic Analysis of a Cube 



Figure 5.14: Cube: geometry and boundary conditions 


The problem of a cube subjected to uniaxial step load and undergoing large deforma- 
tion has been selected for single body elasto-plastic analysis. The geometry and boundary 
conditions are shown in Figure (5.14). The geometric and material properties are: 

• Side of the cube (L) = 30mm 

• Young’s modulous (E) = 1114.58 N/mm 2 

• Poisson’s ratio (i/) =0.36 

• Density (p) = 2000 x 10 -09 Kg /mm 3 

• Hardening coefficient (K)= 87.42 N/mm 2 

• Hadening exponent (n) = 0.909 

• Initial yield stress (°<T y ) = 44.58 N/mm 2 




The discretization used is: 


• Number of elements = 1 (1 x 1 x 1) 

• Number of degree of freedom = 24 

The convergence tolerence used is 1%. .The equivalent stress is plotted against equivalent 
plastic strain in Figure (5.15) for both the Jaumann stress rate and new incremental 
objective stress measures. The actual hardening law 

= % + K (eg,)" (5.1) 

is also plotted in the same figure. The results show that the responses predicted by the 
two different stress measures are the same and match well with the actual hardening curve 
of the material. The unloading is also predicted correctly. 



Figure 5.15: Elasto-plastic response of the cube 




5.3.2 Dynamic Elasto-plastic Response of a Fixed Square Plate 




Figure 5.16: Plate: geometry and boundary conditions 


The impact response of a square plate fixed at all its edges and charaterised by the 
elasto-plastic behaviour is obtained. Figure (5.16) shows the geometry and boundary 
conditions of the plate. Only one quarter of the plate has been modeled due to symmetry 
along X and Y axes. All specified quantities refer to the quarter plate model. The 
geometric and material data for the plate are: 

• Side of the plate (L) = 125mm 

• Thickness of the plate (h) = 25mm 

• Young’s modulous (E) = 1114.58 N/mm 2 

• Poisson’s ratio (i/) = 0.36 

• Hardening coefficient (K)= 87.42 N/mm 2 

• Hadening exponent (n) = 0.909 

• Initial yield stress (°cry)= 44.58 N/mm 2 

• Density of the plate material (p) = 2000 x 10 -09 Kg /mm 3 




Mesh and Time — step Convergence Study 

As stated earlier, mesh and time-step convergence studies are carried out to select an 
appropriate mesh and time step. The mesh and time steps used are given below. 

1. Mesh Convergence: The time step used is At = 0.3 xlO -03 sec and the convergence 
tolerance used is 1%. The discretizations used are: 

• Discretization 1: 

— Number of elements = 200 (10 x 10 x 2) 

— Number of degrees of freedom = 1089 

• Discretization 2: 

— Number of elements = 1600 (20 x 20 x 4) 

- Number of degrees of freedom = 6615 

From the mesh convergence study, it is found that the second discretization gives 
better results. So it is chosen for the final comparison. 

2. Time Convergence: The time steps used are stated below. The mesh used is 
20 x 20 x 4 and the convergence tolerance used is 1%. 

• Time-step 1: 0.6 x 10“ 03 sec 

• Time-step 2: 0.3 x 10 -O3 sec 

It is found that both the time steps give nearly the same results. For the final study, 
both the time steps are used. 

The results obtained by using both the Jaumann stress rate and new incremental ob- 
jective stress measures are shown in Figures (5.17-5.18). Use of Jaumann stress rate pre- 
dicts a lower displacement response (Figure 5.17) but higher stress response (Figure 5.18) 
compared to that of the new incremental objective stress measure. 

This can be explained as follows. It is well-known that the Jaumann stress rate is 
objective only when the rotation of the particle is small. When the particle rotation 
is large, it does not completely eliminate the change in Cauchy stress components due 
to rotation while calculating the incremental stress from the incremental strain. In the 
present problem, the particle rotation is large in most of the domain, especially at the point 
of maximum equivalent stress. This may be the reason for the higher prediction of the 



maximum equivalent stress when the Jaumann stress rate measure is used. However in the 
previous problem (Section 5.3.1), there is no rotation, therefore, there is no discrepancey 
in the two responses. 

To analyse the situation further, the variation of equivalent stress with equivalent plas- 
tic strain is plotted at a typical Gauss point (centre of the plate) for both the discretiza- 
tions (Figures 5.19-5.20). The figures show that, even for a coarse mesh, the predictions of 
the new incremental objective stress measure match very well with the actual hardening 
curve of the material. On the other hand, the predictions of the Jaumann stress rate 
measure match the actual hardening curve only for the fine mesh. For a coarse mesh, 
the level of maximum equivalent stress predicted by the Jaumann stress rate measure is 
higher. 
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Figure 5.17: Maximum displacement response of the plate 



Figure 5.18: Maximum equivalent stress response of the plate 





Figure 5.19: Variation of equivalent stress with equivalent plastic strain (5x5x1) 



Figure 5.20: Variation of equivalent stress with equivalent plastic strain (20 x 20 x 4) 





5.4 Large Deformation Elasto-Plastic Contact-Impact 
Analysis 

5.4.1 Elasto-plastic Impact of a Sphere on a Plate 



Figure 5.21: Impact of a sphere on plate 


The problem of an impact of a sphere with a plate is selected for elasto-plastic contact- 
impact analysis (See Figure (5.21)). The sphere is given an initial velocity, V = (10,0,-10), 
in the coordinate system shown. The plate is fixed along the edges. The coefficient of 
friction used is /i = 0.2. 

The geometric details of the sphere and plate are shown in Figure (5.21). The material 
properties assumed for both the bodies are as follows: 

• Young’s modulus ( E ) = 4 x 10 5 units 

• Poisson’s ratio (v) = 0.3 

• Density ( p ) = 4 units 

• Initial yield stress ((cj,)o) = 5000 units 



• Hardening coefficient (K) = 40000 uni ts 

• Hardening exponent (n) = 0.4 

Time step used for the analysis is 0.002 units. The mesh details of the plate and sphere 
are as follows: 

• Number of elements for sphere = 672 (See Figures (5.22) and (5.23)) 

• Number of elements for plate = 1024 (16 x 16 x 4) 

• Number of d.o.f. for sphere = 2319 

• Number of d.o.f. for plate = 4335 



Figure 5.22: Finite element plot of sphere (Full view) 

The X and Z displacements of the bottom node of the sphere are shown in Figures 
(5.24) and (5.25) for the frictionless case and in Figures (5.26) and (5.27) for the case 
of friction. It is observed that the displacement prediction of both the stress measures 
match well. Convergence difficulties were encountered for the case of friction for the new 
incremental objective stress measure after certain increment number. 

The variation of maximum equivalent stress in both the sphere and plate for the 
frictionless case is shown in Figure (5.28). It is observed that, as in the case of single body 




Figure 5.23: Finite element plot of sphere (Sectional view in X-Z plane) 

analysis, the use of Jaumann stress rate measure predicts the higher value of maximum 
equivalent stress for both sphere and plate. The results for the frictional case are shown in 
Figure (5.29). Here also, the maximum equivalent stress predicted by the use of jaumann 
stress rate is higher. 

The analysis was repeated for sphere and ball made of polycarbonate material with 
similar geometry and mesh details. The material properties of this material are given 
below: 

• Young’s modulus ( E ) = 1.11458 x 10 9 N/m 2 

• Poisson’s ratio (v) = 0.36 

• Density (p) = 2000 kg/m 3 

• Initial yield stress = 4.458 x 10 7 N/m 2 

• Hardening coefficient ( K ) = 8.742 x 10 7 N/m 2 

• Hardening exponent (n) = 0.909 

The velocity given to the sphere is (10,0,-10) m/s, in the coordinate system shown in 
Figure (5.21) and the coefficient of friction at the contact interface is 0.3, for the case of 
analysis with friction. The time step used in this analysis is 0.0003 seconds. The X and 




Z displacements of the bottom node of the sphere for the frictionless case are shown in 
Figures (5.30) and (5.31). The results for the frictional cases are shown in Figures (5.32) 
and (5.33). The variation of maximum equivalent stress in both the sphere and plate is 
shown in Figure (5.34) for the frictionless case and in Figure (5.35) for the case of friction. 
Convergence difficulties were encountered in the frictional case, and therefore the analysis 
was done upto 0.0123 seconds for the case of Jaumann stress rate and 0.0081 seconds for 
the case of new incremental objective stress. As in the earlier cases, the use of Jaumann 
stress rate predicts the higher value of maximum equivalent stress. 



Figure 5.24: Displacement of the sphere along X-direction for frictionless case. 











Figure 5.27: Displacement of the sphere along ^-direction for frictional case. 



Figure 5.28: Ma ximum equivalent stress in sphere and plate for frictionless case 




Figure 5.29: Maximum equivalent stress in sphere and plate for frictional case 



Figure 5.30: Displacement of the sphere along ^-direction (polycarbonate material) fric- 


tionless case 





Figure 5.31: Displacement of the sphere along Z-direction (polycarbonate material) fric- 
tionless case 



Figure 5.32: Displacement of the sphere along X-direction (polycarbonate material) fric- 


tional case 






Figure 5.33: Displacement of the sphere along Z-direction (polycarbonate material) fric- 
tional case 



Figure 5.34: Maximu m equivalent stress in sphere and plate (polycarbonate material) 


frictionless case 








Chapter 6 


Conclusions and Scope for Future 
Work 

6.1 Conclusions 

The objective of the present work is to compare the realative performance of two objec- 
tive stress measures namely the incremental Jaumann stress rate measure and the new 
incremental objective stress measure. The comparison is carried out for the following four 
cases: 

1. Large deformation, dynamic, elastic, single body analysis. 

2. Large deformation, elastic, contact-impact analysis. 

3. Large deformation, dynamic, elasto-plastic, single body analysis. 

4. Large deformation, elasto-plastic, contact-impact analysis. 

The following conclusions can be drawn: 

1. For single-body elastic analysis, both the stress measures give the identical results. 

2. For elastic contact analysis, the maximum equivalent stress predicted by the two 
stress measures is somewhat different. 

3. For single-body elasto-plastic analysis, the use of Jaumann stress rate measure pre- 
dicts higher m aximum equivalent stress compared to the new incremental objective 
stress measure. However, in problems involving no rotation, both the responses are 
identical. 


4. For the case of elasto-plastic contact analysis also, the use of Jaumann stress rate 
predicts higher maximum equivalent stress compared to the new incremental objec- 
tive stress measure. 


6.2 Scope for Future Work 

The desired future extensions of this work are the following: 

1. Currently the programme has the provsion of Jaumann stress rate measure and 
new objective stress measure. Other stress measures along with appropriate strain 
measures can be incorporated and a comparitive study can be performed. 

2. Advanced solution procedures such as the arc length method may be employed to 
handle snap-through and snap-back responses in structures, since such cases cannot 
be handeled by the Newton-Raphson technique. 

3. Finite elements other than eight noded brick element can be incorporated. 

4. Material damping effects can be included. 

5. Option for using penalty function method for developing the contact stiffness matrix 
can be incorporated. 

6. Other material models like kinematic and anisotropic hardening and hardening laws 
other than power law can be implemented. 

7. Non-classical friction laws can be used for modelling the effect of friction. 

8. Finally, it is highly desirable to have versatile pre and post processors, for the 
generation and interpretation of results in 3-D analysis. 
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